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ABSTRACT 


Standard theoretical methods of analysis which work well for seem- 
ingly more complex problems fail to predict the experimentally observed 
instability of fully developed, incompressible pipe flow at any Reynolds 
number. Past research by Harrison and Arnold on the stability of pipe 
flow yielded erroneous results due to errors in the setup of the problem 
and formulation of the boundary conditions at the axis. 

A revised theory with particular attention to the rather complex 
boundary conditions at the axis has recently been developed. Improved 
finite differencing techniques with consistent fourth order truncation 
error were also used to approximate the governing differential equations. 

Numerical resuits for angular wave numbers zero and six show that 
the flow is stable at al! Reynolds numbers. Results for angular wave 
number one contain instabilities at al] Reynolds numbers for small 
values of the axial wave number. These results are tabulated, plotted, 
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I. INTRODUCTION 


Osborne Reynolds [Ref. 1] conducted his classical experiments on the 
transition from laminar to turbulent flow of fluids in circular pipes 
nearly 100 years ago. Based on his experiments the critical Reynolds 
number for pipe flow was found to be 2300. An analytical solution for 
the problem of predicting instabilities of fully developed, three dimen- 
sional, incompressible flow of constant viscosity fluids in pipes has 
been pursued actively since then. Several approaches have been taken in 
solving the inherently nonlinear Navier-Stokes equations, which along 
with the continuity equation, govern the behavior of fluid flow in pipes 
of circular cross-section. 

Salwen and Grosch [Ref. 2] studied pipe flow with purely sinusoidal 
streamwise (axial) perturbations. Infinitesimal velocity and pressure 
perturbations, which were explicit functions of time and not complex or 
exponential in form, were introduced into the flow field. Numerical 
calculations were carried out at various angular wave numbers and it was 
concluded that the flow was stable for all] axial wave numbers and 
Reynolds numbers. Perturbations with exponential growth in space but a 
purely sinusoidal variation in time were explored by Garg and Rouleau 
[Ref. 3] and found to be stable. Gill [Ref. 4] studied combinations of 
exponential growth in space and in time using power series analysis and 


concluded that these flows were all stable. 
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Because of the general agreement that pipe flow is stable to infini- 
tesimal disturbances, two other approaches to the problem of predicting 
instabilities at the experimentally observed critical Reynolds number 
have been pursued. First, investigations by Davy and Drazin [Ref. 5] 
confirmed that pipe Poiseuille flow was stable at all Reynolds numbers 
for infinitesimal disturbances that were both temporally and axially 
complex and exponential in nature. However, despite the fact that the 
governing equations must necessarily remain nonlinear with the introduc- 
tion of disturbances of finite amplitude, they concluded that the flow 
is unstable with respect to finite disturbances. Similar theoretical 
examinations of plane Poiseuille flow between two parallel plates led 
McIntire and Lin [Ref. 6] to this same conclusion. Second, it was 
postulated by Huang and Chen [Ref. 7] and Leite [Ref. 8] that the origin 
of flow instabilities occur in the entrance region of the pipe, where 
the flow has not yet become fully developed. Both theoretical and 
experimental studies have been done to support these hypotheses. Garg 
{Ref. 9] also found that flow instabilities existed near the entry 
region for developing pipe flow. He concluded that results using the 
Hornbeck velocity profile more nearly approximated experimentally deter- 
mined values of the critical Reynolds number. While these investi- 
gations have shown instabilities to exist in pipe flow, a completely 
generalized solution to the linearized problem of fully developed pipe 
flow has never been accomplished. 

Recently, a more general theory consisting of perturbations which 


have a fully complex exponential form with respect to time and the axial 
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coordinate, a purely imaginary expenential form in the angular coordi- 
nate and are three dimensional in nature have been explored. Develop- 
ment of this theory by Harrison [Ref. 10] and further numerical in- 
vestigations by Arnold [Ref. 11] failed to produce conclusive results 
that instabilities exist in fully developed pipe flow. Errors in the 
setup of the problem and inadequate formulation of the boundary condi- 
tions at the pipe axis contributed to these erroneous results. Arnold, 
however, was on the right track in his research by incorporating newly 
formulated boundary conditions at the axis developed by Gawain [Ref. 12]. 

Arnold's work, for the most part, remains valid except that he 
overlooked a small but significant detail that resulted in unclear 
definitions of stability and instability. He used a complex axial wave 
number a, of the form a = Ap * jas, whereas the correct version should 
simply be a purely imaginary quantity of the form ia. Thus; Arnold 
introduced an extra degree of freedom, namely the fictitious and in- 
correct quantity Op. This accounts for the erroneous instabilities 
which he obtained in his numerical investigation of pipe flow. 

Further work by Gawain [Ref. 13] refined the boundary conditions at 
the axis and incorporated the previously defined imaginary axial wave 
number, in exponential form, which represents only sinusoidal osciila- 
tions with respect to the axial coordinate. Using these advancements in 
the linearized theory, Gawain showed an asymptotic trend toward neutral 
stability of the flow with increasing Reynolds numbers for angular wave 
number, n = 0. This paper duplicates and confirms the results for 
n= 0, that the flow is stable at all Reynolds numbers and axial wave 


numbers. It is further shown that fully developed pipe flow for n = 6 
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is also stable and shows the same trends as for n = 0. Instabilities 
are shown to exist in the numerical analysis of the vorticity transport 
equations for n = 1. The paradox that now appears is that for angular 
wave numbers n = 0 and 6 the flow is apparently stable at all] Reynolds 
numbers while for n = 1 the flow is apparently unstable at all Reynolds 
numbers. Neither cf these theoretical results is consistent with the 
known experimental fact that there exists a critical Reynolds number 
below which pipe flow is stable and above which it is unstable. However 
the new result for n = 1, while it can hardly be called correct, is at 
least encouraging in that it shows for the first time that instabilities 
can in fact exist in the solution of the linearized vorticity transport 
equations for fully developed pipe flow. 

In addition to implementing the improved axis boundary conditions 
and the advancement in the linearized theory, improved finite differ- 
encing techniques were used to approximate the vorticity transport 
equations along a uniform radial mesh. Highly accurate numerical solu- 
tions for the previously mentioned angular wave numbers, n = QO, l, 
and 6, were obtained using double precision, complex numbers. Addition- 
ally, finite difference equations with consistent fourth order trunc- 
ation error were used throughout. Details of the development of the 
vorticity transport equations and the central and non-central finite 
difference equations are discussed in later chapters. Systematic and 
extensive calculations remain to be accomplished for anguiar wave 
numbers, n = 2, 3, 4, 5, and n > 6. The results of these numerical 


analyses should prove to be fruitful in establishing a theory which 
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adequately explains the experimentally observed fact of the onset of 


flow instabilities at the critical Reynolds number. 
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, II. THE VORTICITY TRANSPORT EQUATION 


A. BASIC THEORY 
The governing equations for laminar flow of fluids in a circular 


pipe are derived from Newton's second law of motion, 
F=ma (2-1) 


and the conservation of mass principle through an infinitesimal control 
volume. Equation (2-1) results in the familiar Navier-Stokes equation, 
which is a vector equation in three cylindrical coordinates and in time. 
The Byimition for conservation of mass results in the continuity equation 


for incompressible flow and is given by, 
"7-V=9 (2-2) 


For the case at hand, the following assumptions are made, 
1. The fluid is incompressible, p = constant 
2. The fluid has constant viscosity, uw = constant 
3. The flow is Tully developed, U = f(r) only 
4. The flow is three dimensional in nature 
5. The mean flow is steady, aV/at = 0 
6. The perturbation flow is unsteady, av/at # 0 
7. The effects of body forces are negligible 
The resulting Navier-Stokes equation can be stated as the sum of the 


forces per unit mass is equal to the acceleration or, 
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The forces acting on the fluid are pressure forces and viscous forces. 
Body forces and hydrostatic forces balance out and are not a factor 
here. 

Equation (2-3) in its dimensional form can be expressed as, 

- yr se 4: (He) y2x Vx = dVv* (2-4) 
oer es dt 

where the starred quantities are fully dimensional. 

It should be noted here that the equations presented through equa- 
tion (2-18) represent the mean flow. Equation (2-4) can also be 


expressed in a nondimensional form as, 


e ee ee = 
Rate VON = Ge (2-3) 
where Re is the Reynolds number and is defined as, 
RU 
BS rE (2-6) 


This definition of Reynolds number is based on the pipe radius R as the 
characteristic length, the mean volumetric velocity U as the character- 
istic velocity, and the kinematic viscosity v which is constant. The 
three quantities on the right side of equation (2-6) are dimensional, 
but Re is a dimensionless quantity. On this basis, the critical 
Reynolds number for transition from laminar to turbulent flow becomes 


1150 (vice the value of 2300 which is obtained if the Reynolds number is 
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defined in terms of pipe diameter). All subsequent equations presented 
in this paper are in nondimensional form. 
Returning to the derivation of the mean flow, the acceieration term 


or substantial derivative, dV/dt can be expanded to the form, 

Car?) 

but (V ° v)V can be further expanded by a well known vector identity to, 
@-wW=vh)-tx xh (2-8) 

Additionally, v2V from equation (2-5) can be expanded to, 


vV=av(v-V)-Vx (Vx V) (2-9) 


ay. - Vv = 0, since continuity must be satisfied and equaton (2-9) 


becomes, 


— 


v2V=-Vx (Vx V) (2-10) 


With appropriate substitutions, the final form of the Navier-Stokes 


equation becomes, 


<yR.-< 2 Wa7 2-3 Can 4s ot : 
VP Re Vx (VV x ¥)= 9 G er SEY On oe 5t (2-11) 
A more convenient form of equation (2-11) is, 
4 ae pa ye A > av y 
V (P + 3 y) = Ra VX AV XV) = VX YK VD + St (2-12) 
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As previously discussed, all the quantities in the equations nave 
been nondimensionalized. Dimensionless cylindrical coordinates x, r, 


and @ are used with e é. and & denoting the unit vectors in the three 


x? 
coordinate directions, respectively. It is well known that the velocity 
profile for fully developed, three dimensional, incompressible flow in 


circular pipes can be expressed analytically by, 
y = 2. (1) = re) (2-13) 


This function has the shape of a paraboloid of revolution as shown in 


Figure 2-1. 


Figure 2-1. Velocity Profile of Fiuid Flow in a Pipe 


The nondimensional mean velocity vector V(r) of the flow is given by, 


Vir) =Ue =2(1-r2)e (2-14) 
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A second vector quantity 6, the mean vorticity, is defined as the cur] 


of the mean velocity vector and is given by, 


which will be used later. 
B. THE VORTICTY TRANSPORT EQUATION 


The vorticity transport equation is now obtained by taking the curl 


of the Navier-Stokes equation (2-12), resulting in, 
2 x = 
~Vxv (P+ (S)) = Vx (Vx (VX H)) - Vx Wx (Vx) + oe V) 


(2-16) 


The primary reason for taking the curl of the Navier-Stokes equation © 
is to eliminate the unknown scalar on the left side of equation (2-16). 
Since the curl of the gradient of a scalar is equal to zero, equation 


(2-6) reduces to, 


Tx (Vx (Vx V))-Vx Wx (Wx ¥)) + eee =9 (2-17) 


i 
Re 


Equation (2-17) can be simplified by substituting the mean vorticity 


vector relation equation (2-15) into equation (2-17) resulting in, 


~ 7x x 6) = vx (¥ x 4) + 2 =0 (2-18) 
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In order to analyze this flow field and in particular to determine 
the transition from laminar to turbulent flow, a disturbance of small 
amplitude is introduced into the flow. It is assumed that the resulting 
flow is made up of the steady state, laminar flow and a small perturba- 
tion flow superimposed on the laminar flow. The total velocity and 


vorticity respectively become, 


Wate eign 
and 
G°-=8+0 (2-20) 


where v and w are the velocity perturbations and vorticity perturbations, 
respectively. It is of course a necessary requirement that the contin- 


uity equation for the total flow be satisfied here as well, 
Vie VoumiWien Wi 2-36) Vor vis 0 (2-21) 


In view of this reiation and of equation (2-2), the perturbation flow 


independently satisfies the continuity condition and 
Ve va= 20 (2-22) 


The introduction of the perturbation quantities into the vorticity 
transport equation is accomplished by substituting V’ for V and 4° for @ 


into equation (2-18). This results in, 


1 


me VX (VX G+w))-Vx(V+tv)x @+u)] + ee 2 


0 
(2-23) 
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which is the equation for the total or resultant flow field. Upon ex- 
panding equation (2-23), applying the mean steady flow assumption, 
aV/at = 0, and subtracting equation (2-18) from it, equation (2-23) 


becomes, 


& Vx (Vxw)-Vx (Wx) - @xv) + Wx w)] + = bl 


(2-24) 


Now equation (2-23) is linear in the perturbation quantities except for 
the second order term, (Vv x w). Since the perturbations are assumed to 
be small, this term can be neglected. 


The perturbation vorticity is defined as, 
Lo aa (2-25) 


which is analogous to equation (2-15) for the mean flow. Making this 
substitution in equation (2-24) yields the linearized vorticity *rans- 


port equation, 


BV X (VK (VX Y)) - 0x Kx) HTK Tx KY) + 0x Be9 
(2-26) 


As previously mentioned, the continuity equation (2-22) for the per- 
turbation velocity must be satisfied. A convenient way to insure this 


is to express vV in terms of a velocity vector potential function W as 


follows, 


; 7 
| p : =< : 
<<g coq! fer wor? sad tuzon, 7 j = 
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v=vxw (2-27) 


It can be verified that if v is defined in this way, equation (2-22) is 
satisfied identically for any arbitrary vector function W. 

Instead of proceeding directly with the real vector functions Vv and 
W, it is advantageous to work with the Fourier transforms of Vv and W, 
which are in general complex functions. This is the form which is used 
here. It can be shown, because of the linearity of the basic equation, 
that if a solution can be found for the Fourier transforms, then the 
corresponding real vector functions will satisfy the basic equation as 
well. A detailed explanation of these functions was presented more 
elegantly by Gawain [Ref. 13]. The vector potential function W, now 


representing the Fourier transform, can be expressed as, 
W=([F(r) é€, + G(r) e+ H(r) é,] e° (2-28) 


where X is a convenient abbreviation for the complex quantity, 
X = lax + ind + yt (2-29) 


The axial wave number a is a real quantity greater than or equal to zero 
which in the exponential form represents purely sinusoidal axial vari- 
ations of v. The angular wave number n is a positive integer and repre- 
sents sinusoidal angular variations of Vv which are periodic with respect 
to coordinate 9. The complex frequency y contained in the exponent is 


of the form, 
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The quantity Yp clearly represents the exponential time rate of growth 
or decay of the perturbation amplitude. The basic criteria for hydro- 
dynamic stability can be defined in terms of YR: Positive values of Yp 
represent growth and signify flow instability while negative values of 
Yp represent decay and signify flow stability. 
The vorticity transport equation can now be expressed in terms of 
> 


the velocity vector, V and the complex vector potential, W. By making 


the substitution for v, equation (2-26) now becomes, 


a VX (VX (Vx (Vx H))) - 0x HK (Wx (x H))) 


+Vx (Vx (Vx (Vx W))) +x (Vx a) =( (2-31) 


In convenient shorthand notation, equation (2-31) becomes, 


Equation (2-31) = F(r)e* = [T.(r) e. + F_(r) e, +T,(r) e ] eX = 0 
x r r 8 8 
(2>32) 
and can be expressed in matrix notation as, 
r 0 
x 
i= ue =A 0 (2-33) 
Tp 0 
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Therefore, the three components of the vorticity transport equation are 


equivalent to three simultaneous scalar equations, namely, 


ie = 0 
ee = 0 (2-34) 
U5 = 0 


Of these three equations, only two are actually independent since, 
as can be seen from equation (2-31), f e% is itself the curl of another 
vector. Therefore f e*~ must necessarily be nondivergent. It can be 


shown that the components of i satisfy the identity, V - (fT e*) SO OF, 


, ut in = = 
lol at r D (rr) + =) U9 = 0 (2-35) 


where D is the symbol for the differential operator d/dr. Two linearly 
independent equations can now be formulated from equations (2-34) given 
the constraint. of equation (2-35). This can be reduced to some appro- 
priate linear combination of two of the equations which is set equal to 
zero and to the third equation which is also set equal to zero. This 
results in a system of two coupled, fourth order differential equations 
in three unknowns. Recall that the vector potential function W is 
expressed in terms of F(r), G(r), and H(r). It can be shown that one of 
these variables is redundant and may be eliminated without any loss of 
generality in the solution of the equations. For the present case, the 


solution assumes its most convenient form if, 


P(r) = 0 (2-36) 
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(9) 


This convention has therefore been adopted and now G(r) and H(r) become 
the two coupled eigenfunctions of the governing fourth order differen- 
tial equations. Returning to the problem of selecting the proper linear 
combination of equations comprising ise it is important to choose this 
combination in an optimum fashion by analyzing the algebraic structure 
of these equations. Table 2-1 is taken from reference [13], corrected, 
and duplicated here for that purpose. Notice that equations (1) and (3) 
in the table are third order in G. If the particular linear combination 
of these equations is chosen as appears in equation (4), the result 
reduces to an equation of second order in G. This does not affect the 


order of the result since the function H remains fourth order. 


Table 2-1. Properties of the Vorticity Transport Equations 


HIGHEST 
TERM OF TERM OF NEGATIVE 
HIGHEST HIGHEST POWER 
EQUATION ORDER IN G ORDER IN H he 
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The formulation of the problem is simplified greatly by choosing two 


independent vorticity transport equations in the form, 


rr =0 (2°37) 


“(Ge f * r, = 0 (2-38) 


C. THE VORTICITY TRANSPORT MATRIX EQUATION 


While the fourth order linear vector operations and the associated 
algebra becomes very lengthy and rather tedious, the investigator at- 
tempting to expand the final form of the linearized vorticity transport 
equation (2-31) should pay sufficient attention to these details. 
Recall that the curl operation on a vector in cylindrical coordinates 


takes the form, 


i 11> > 
r ey r SP 6 

yxW= | & = ~ (2-39) 
0 Ger rHe* 


After all the nabla (V) operations are complete, the final equation wil] 
be in the form of a vector equation as defined in equation (2-32). The 
first vorticity transport equation consists of the coefficient repre- 


sented by r .. and is set equal to zero. The second equation consists of 
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the previously defined linear combination of the coefficients represent- 
ed by cS and lg and is also set equal to zero. The two resulting 
simultaneous ordinary differential equations are then expanded and the 
coefficients of the derivatives regrouped so that equations (2-37) and 


(2-38) can be conveniently expressed in matrix format as follows, 


D*G (03a) D2G 0G ) 
ma + [(M,] eal + [M,] > ST | 


D4H DH oH | 
G D2G {0G (¢) 

+ [Mo] | = s(t» + [Ny] | eA 3 ) (2-40) 
H D2H l oH lH | 


As it turns out, the equations are a pair of coupled homogeneous differ- 
ential equations that can be solved as an eigenvalue problem. This 
technique is discussed in the chapter on Numerical Methods. 

The matrices which appear in equation (2-40) are 2 X 2 matrices and 
are summarized in detail in Appendix A. Equation (2-40) and the re- 
spective matrix elements contained in Appendix A are the generalized 
vorticity transport equations for all] angular wave numbers n, at all 
axial wave numbers a, and all Reynolds numbers. As will be explained in 
detail in the chapter on Boundary Conditions, the generalized equations 
are transformed by appropriate changes of variable depending on the 
angular wave number. The changes of variables from the eigenfunctions 
G(r) and H(r) to the new eigenfunctions P(r) and Q(r), respectively, are 
determined by the original boundary conditions that must be satisfied, 


particularly at the pipe axis. 
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D. EQUATIONS FOR n= 0 


For the case n = 0, equation (2-40) is greatly simplified in that 
the two equations uncouple and allow an independent investigation of 
either eigenfunction, G(r) or H(r). Miscellaneous trial calculations 
made earlier suggest that the solutions for G(r) are probably stable for 
all values of @ and Re. _ Therefore, only the solutions for H(r) were 
investigated. It should also be noted that the examination of the axial 
perturbation velocity is of particular interest in this analysis and, in 
the case for n = 0, is not a function of G(r). The change of variable 


required to satisfy the boundary conditions is, 
HCr ) f= er sQGr) (2-41) 
Taking the derivatives of H(r) yields, 


DH 


r DQ + Q 


D2H =>r D*Qc+ 2.0Q 


D7H =r D2Q + 3: D2Q 


D*H = r D4Q + 4 D3Q (2-42) 


With the substitution of equations (2-41) and (2-42) into equation 


(2-40), the following expression is obtained, 


Ma DQ + Mg D3Q + Mo D2Q + Mz DQ + Mo Q = y (Ng D2Q + Ny DQ + No Q) 


(2-43) 
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The primed coefficients are actually the elements from row 2 and column 
2 of the newly transformed coefficient matrices and are not enclosed in 
matrix brackets. However, in general, the new coefficient matrices [M-] 
and (NGI, (where the indices 1 = 0, 1, 2, 3, 4 and j = 0, 1, 2), are 
formed after the change of variable is made and the terms are coi lected. 
The primed coefficients used to solve the transformed vorticity trans- 
port equation in terms of Q(r) for n = 0 are defined in Appendix B, Part 


A. 


E. EQUATIONS for n 


i} 
r=) 


For the case n = 1, eigenfunctions G(r) and H(r) do not become un- 
coupled and must be solved for in a system of simultaneous equations in 
a form similar to equation (2-40). An additional parameter H(0) is 
introduced into the equations by the change of variables. This parti- 
cular change of variables is required because of the complicated nature 
of the boundary conditions at the axis and is thoroughly explained in 
the chapter on Boundary Conditions. For n= 1, the change of variables 


required to satisfy the boundary conditions are, 


G(r) = -7. HOO) + 72° PCr) 
(2-44) 
HGr = HCO) + 72 00r) 
Taking the derivatives of G(r) yields, 
DG = r2 DP + 2r P 
(2-45) 


D?G = r2 D2P + 4r DP + 2 P 
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and taking the derivatives of H(r) yields, 


DH oS 72 DQ + 2rQ 
Deny= £2,0°Q + 4r 00 + 2 Q 

(2-46) 
Be = r* D7Q + Gr 070 + G6 OQ 


04H = r@ D*Q + Sr DQ + 12 02Q 


Remember that in general D(H(0)) # DH(O), where DH(O) is the value of 
the first derivative of H(r) evaluated at the point r = 0. Since 
D(H(O)) = 0 here, it vanishes from the expressions for the derivatives 
of G(r) and H(r). In order to accommodate H(0) in the system of equa- 
tions where it appears explicitly, two additional 2 X 1 column matrices 
are required, namely [Ms] and [Nj]. After the changes of variable have 


been made, equation (2-40) now becomes, 


ee OP (PPP fF 
a. i gs bres al Hist gine’ tal 


D2Pp oe P } 
+ [MZ] {H(0)} = y («ws | ede <)> Ned) te TNS 1H(0)1) 
n2qJ (naj lo} 


(2-47) 


The coefficient matrices of equation (2-47) used to solve the trans- 
formed vorticity transport equations in terms of P(r) and Q(r) for n= 1 


are defined in Appendix B, Part B. 
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F. EQUATIONS FOR n = 6 


For the case n = 6, eigenfunctions G(r) and H(r) remain coupled. 
The change of variables does not introduce any additional parameters and 
the system of simultaneous equations is solved in the form of equation 
(2-40) once again. The change of variables required to satisfy the 
boundary conditions for n = 6 and additionally for all angular wave 


numbers n > 6 are, 


G(r) 
H(r) 


re P(r) 


r= Q(r) (2-48) 


Taking the derivatives of G(r) yields, 


0G r= DP°+-4r> P 


D2G = r* D2P + Sr? DP + 12r2 P (2°29) 


and taking the derivatives of H(r) yields, 


DH =r? DQ + 3r2 Q 
D2H = r? D2Q + 6r2 DQ + Gr Q 
D?H = r? D3Q + 9r2 D2Q + 18r DQ + 6 Q 


D4H = r® 07Q + 12r? D3Q + 36r 02Q + 24 DQ (2-50) 


Substituting equations (2-48), (2-49), and (2-50) into equation (2-40) 
results in the transformed coefficient matrices for the system of 
equations. The coefficient matrices of the functions P(r) and Q(r) and 
their respective derivatives for n = 6 are defined in Appendix B, 


Part C. 
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III. BOUNDARY CONDITIONS 


The incorrect formulation of the boundary conditions for the linear- 
ized vorticity transport equations at the pipe axis has probably been 
the primary reason that previous investigators have failed to predict 
the onset of flow instabilities in circular pipes at any Reynolds number. 
Gawain [Ref. 13] suggested that if these boundary conditions are formu- 
lated in a rigorous and systematic fashion and applied to the vorticity 
transport equations, a correct numerical solution may be computed which 
agrees with experimental results. It can be seen from equations (2-40) 
and from the coefficient matrices of Appendix A that the two vorticity 
transport equations are coupled, except for the case n = 0. The re- 
sulting pair of differential equations are second order in G(r) and 
fourth order in H(r). In order to obtain a determinate solution of the 
equations, a total of six boundary conditions are required. More 
specificaliy, at least two of the boundary conditions must involve G(r) 
or its derivatives and at least four of the boundary conditions must 
involve H(r) or its derivatives. It might be expected that when the 
equations are coupled, the boundary conditions may also be coupled. 
That is to say that the boundary conditions may consist of various 
linear combinations of the functions and their derivatives, which turns 


out to be the case and is shown later. 
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The boundary conditions at the pipe wall are determined in a fairly 
straight forward manner. They are derived from the fact that the com- 
ponents of the perturbation velocity must vanish at the wall. From the 


definition of the perturbation velocity, defined previously as, 
v=u x (2-27) 


it can be shown that the following three boundary conditions result. 


G{1) = 0 
Hi), = .0 
DH(1) = 0 (3-7) 


Since six boundary conditions are required and three have been deter- 
mined at the pipe wall, there remain three boundary conditions to be 
determined at the pipe axis. Of these three conditions, at least one 
must involve G(r) or its derivatives and at least two must involve H(r) 
or its derivatives. 

The proper derivation of the three boundary conditions at the axis 
turns out to be a non-trival problem. As can be seen in Table 2-1, the 
highest negative power of r that appears in the equations is ie 
Closer inspection of the vorticity transport equations of equation 
(2-40) and the coefficient matrices in Appendix A reveals that the 
equations contain terms in ae ot rsiea ana and NG At first glance, 
it may appear that these equations are not satisfied in the limit as r 
approaches zero. It should be noted that it would be incorrect to 


multiply the equations through by any positive power of r since that 


would alter the degree of the singuiarity of the equations at the axis. 
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This problem can be resolved and subsequently the boundary conditions at 
the axis can be deduced rigorously by expanding the functions G(r) and 
H(r) as power series in r. The highest negative power of r is fo and 
appears in the coefficient matrix [Mj] for the pair of functions G(r) 
and H(r). The series expansion for G(r), and similarly for H(r), is 


carried to the fourth derivative as follows, 


r2 r3 r4 
G(r) = G(O) + DG(O)r + D2G(0) nat D°G(0) =, + D4*G(0) a ee (3-2) 


It can be seen that when the series expansion expressions for G(r) and 
H({r) are substituted into equation (2-40), the higher order terms of 
equation (3-2) that contain powers of r greater than four will vanish in 
the limit as r approaches zero. For the first derivatives of G(r) and 
H(r), the highest negative power of r that appears in the coefficient 
matrix [M,] is oo Therefore, the series expansion for the first 
derivative of G(r), and similarly for H(r), up to the fourth derivative 


term is, 


DG(r) = 0 + DG(O) + D2G(0)r + D2G(0) c + 046(0) . oe (3-3) 


The remaining derivatives of the functions are expanded in a similar 
manner through the fourth derivative term. The resulting power series 
approximations for G(r) and H(r) and their derivatives are then substi- 
tuted into equation (2-40). The coefficients are then regrouped in 


ascending powers of r. 
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As mentioned previously, the equations must be satisfied in the 
limit as r approaches zero. Therefore, the coefficients of the func- 
tions that contain “ “ans eae Sp and or must all be zero. The 
remaining coefficients will contain positive powers of r and will vanish 
in the limit as r approaches zero. By considering only the terms con- 
taining negative powers of r and setting those coefficients equal to 
zero, five pair of equations in ten unknowns remain so that the boundary 
conditions can be determined. The unknowns are the functions of G(r) 


and H(r) and their first four derivatives evaluated at the pipe axis, 


r=06. It is convenient to summarize these equations in matrix format 


below 

Ba al te - 
Bol “fo 5: 
Svea “fol MEET oh “fal 


n(n2=4) $0266] lal rhe, OC ae 
Re [Cs] losucoy) Zia \1 y} [Ce] { onc0) | lo (3-7) 
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| p4Hc0) | | o2Hc0) { 


Kast be 
+ Cia [Cg] - ya? [Dg]) = (3-8) 
luo) fof 
The 2 X 2 matrices which appear in equations (3-4) through (3-8) are de- 
fined in detail in Appendix C. Special conditions at the axis arise 
when the determinants of these matrices are evaluated at specific 
angular wave numbers. They are summarized in Appendix D. 

For each specific value of angular wave number n, it is possible to 
derive a set of boundary conditions at the axis. Special conditions at 
the axis occur for angular wave numbers that cause the coefficients in 
equations (3-4) through (3-8) to equal zero or the determinants in 
Appendix D to be zero. That is e say, if a coefficient vanishes or the 
pair of equations are linearly dependent, the functions may be arbi- 
trarily specified and are in general not equal to zero. But if the 
coefficient, appearing in equations (3-4) through (3-8) is non-zero and 
the determinant exists, the pair of functions must be identically zero. 
This situation occurs for all five pairs of variables when the angular 
wave number n 2 6. 

By taking these special conditions into account for the case n= 0, 


it can be deduced from equations (3-6) and (3-8) that, 


(G(0) } _ f02Gco)) | | Ae 


lucoy) — [ozHcayf ~~ lo 
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Carrying out this procedure for n = 1, 2, 3, 4, 5, and 6 yields a sim- 
plified set of boundary conditions at the axis. The following sets of 
equations (3-10) through (3-16) are essentially the same as those of 
Gawain [Ref. 13] but were independently checked and appear here with a 


few minor corrections. 


n=0 
G(0) = 0 H(0) = 0 (a) aa 
D2G(0) = 0 D2H(0) = 0 (b) 
n=1 
G(O) + 1 H(0) = 0 
(a) 
DG(0) = 0 DH(O) = 0 
D?G(0) = 0 D3H(0) = 0 
3-11 
: {0*G(0)] (026(0)) 
so [C7] + Cia [Cg] - ¥ (Ds) (b) 
[o#Hc0)| [p2H(0)| 
| jec0)) _ {0} 
= (7a [Co] - ye? [05)) | aco) | = lo 
n=2 
G(0) = 0 H(O0) = 0 s 
a 
DG(O0) + i DH(0) = 0 
a-----------------+---------------------------------------- (3-12) 
02G(0) = 0 D2H(0) = 0 
(b) 
D4G(0) = 0 D4H(0) = 0 
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G(0) = 0 H(0) = 0 
(a) 
DG(0) = 0 DH(O) = 0 
D2G(0) + i D2H(O0) = 0 
(b) 
D2G(90) = 0 D7H(0) = 0 
n=4 
G(0) = 0 H(0O) = 0 
DG(0) = 0 DH(O) = 0 (a) 
D2G(0) = 0 D2H(O0) = 0 
D2G(0) + i D3H(0) = 0 
(b) 
D*#G(0) = 0 D*H(O) = 0 
n=5 
G(0) = 0 H(0) = 0 
DG(0) = 0 DH(O) = 0 (a) 
D2G(0) = 0 , BFR) = 0 
D7G(0) = 0 D3H(O0) = 0 
(b) 
D*G(0) + 7 D*#H(0) = 0 
n2 6 
G(0) = 0 
DG(O0) = 0 H(O) = 0 
D2G(0) = 0 DH(O) = 0 (a) 
D2G(0) = 0 D2H(0) = 0 
D*G(0) = 0 D2H(O0) = 0 
(b) 
D*H(O) = 0 
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Recall that exactly three boundary conditions are required at the 
axis to finalize the solution of the vorticity transport equations. 
Equations (3-10) through (3-16) contain from four to ten boundary con- 
straints including linear combinations of several boundary conditions, 
as suggested earlier. By introducing appropriate changes of variables, 
the pair of functions in G(r) and H(r) and their respective derivatives 
can be transformed to a new pair of functions in terms of P(r) and Q(r) 
and their derivatives. The transformed system of equations take the 
form of equation (2-40) for all values of n, except for n= 1 where 
equation (2-47) is applicable. The change of variables is dependent on 
each specific value of n up through six and is chosen in such a way that 
the boundary conditions in subset (a) of equation (3-10) through (3-16) 
are satisfied identically. The boundary conditions in subset (b) must 
be satisfied explicitly in the finite difference equations near the 
axis. Development of the finite difference equations is discussed in 
the chapter on Numerical Methods. 

The following sets of equations are taken from Gawain [Ref. 13] and 
are included here so that a complete treatment of the subject on bound- 


ary conditions is contained in this paper. 
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Change of Variables 


G(r) = rP(r) 
(a) 
H(r) = rQ(r) 
Boundary Conditions at Axis 3°17} 
DP(O) = 0 DQ(0) = 0 
(b) 
03Q(0) = 0 
Boundary Conditions at Wal] 
P(i)c=n0 Q(1) =0 
(c) 
DQ(1) = 0 


The solutions for P(r) and Q(r) become uncoupled for n = 0 


n=1 


Change of Variables 


G(r) = -1 HCO) + r2 P(r) 
Atr) = HCO) + r2 Q(r) pa 
Sei eaaesat iis 0U[ClC( i(i(i‘“ SC*‘(; (3-18) 
DP(0) = 0 DQ(0) = 
ee ic) 02P(0)) oe fPco)) ae {-i we 
(029¢0)| (ac0)} (a 


Boundary Conditions at Wal] 


“2 A(0) + PC1).= 0 H(O) + Q¢1) 


I 
(2) 


(c) 
-2H(0) + DQ(1) 


i} 
j=) 


ad 


Ch =e 


a 


zi eh 328 onereihnad 


c= (i,t 


wel tuted ae 


n=2 


Change of Variables 


G(r) = - i r DH(O) + r2 P(r) 


H(r) = r DH(O) + r2 Q(r) 


Boundary Conditions at Axis 


P¢o)) = 0 Q(0) = 0 
D2P(0) = 0 D2Q(0) = 0 
Boundary Conditions at Wal] 
=" TADHCO) + PCL) = 0 DH(O) + Q(1) 


~DH(O) + 0Q(1) 


G(r) = r2 P(r) 
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Boundary Conditions at Axis 


P(O) + i Q(0) =0 


Boundary Conditions at Wal] 
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n=4 


Change of Variables 


G(r) = r? P(r) 


(a) 
Ath) = fF O(r) 
Boundary Conditions at Axis (3-27) 
PCO) + i QO) = 6 
(b) 
DP(O) = 0 DQ(0) = 0 
Boundary Conditions at Wal] 
P(1) =0 Q(1) = 0 
(c) 
DQ(1) = 0 
n=5 
Change of Variables 
G(r) = r? P(r) 
(a) 
BT 2:.Qtr) 
Boundary Conditions at Axis (3-22) 
P(O) = 0 Q(o) = 0 
(b) 
DP(O) + i DQ(O) = 0 
Boundary Conditions at Wall 
PCE): =10 Q(1) = 0 
(c) 
DQ<2) = 0 


46 


G(r) = r# P(r) 
(a) 
Cr) = OCP) 
Boundary Conditions at Axis (3-23) 
P(O) = 0 Q(0) = 0 
(b) 
DQ(0) = 0 
Boundary Conditions at Wal] 
P(1) = 0 Q(1) = 0 
(c) 
DQ(1) = 0 


It should be noted that the change of variable relations in subset (a) 
of equations (3-17) through (3-23) does not change the order of the 
derivatives appearing in the transformed vorticity transport equations. 
Upon close examination of subsets (b) and (c) of these equations, it is 
apparent that the correct number of boundary conditions required for the 
solution of the transformed equations has resulted from the change of 
variables. The vorticity transport equations are now expressed in terms 
of P(r) and Q(r) and their derivatives and have exactly three boundary 
conditions at the pipe wall for all values of n and three boundary 
conditions at the axis for all values of n, except for n= 1 and 2. The 
fact that there are four boundary conditions at the axis for angular 


wave numbers n = 1 and 2 warrants further comment here. 


47 


—- 


st vA F5 ars: 


= 7 | 
© yo! simp 


Wier 
if T Ae 
] LO. 
Tsboued 
he” 


‘ a fe Se 1 } 4 th - 
; Lark om u ¢ 


Notice in equations (3-18a) and (3-19a) for the cases n = 1 and 2 
that an additional parameter is introduced by the change of variable 
equations. The additional parameter for n = 1 is H(Q) and for n = 2 is 
DH(O). The introduction of the parameter for each case requires the 
specification of four boundary conditions at the axis instead of the 
usual three. This conclusion was reached by Gawain [Ref. 13] and is 
analyzed in detail in his paper. There is one other peculiarity about 
the boundary conditions at the axis for n = 1 that should be mentioned. 
The axis boundary conditions in equation (3-18b) contain a pair of 
equations that involve the complex perturbation frequency y, which 
represents an eigenvalue of the solution for the equations. The pres- 
ence of y in the axis boundary condition for n = 1 is an important 
factor in formulating a solution for the vorticity transport equations 
and is included in the eigensystem of the finite difference equations. 
The details of the implementation are discussed later. 

The solution of the problem for predicting the onset of flow in- 
stabilities in circular pipes is now at hand. The original vorticity 
transport equations (2-40) are transformed to the functions P(r) and 
Q(r) and their respective derivatives by the appropriate change of 
variables depending on the specific angular wave number being investi- 
gated. After imposing the rigorously deduced boundary conditions, a 


solution can be determined in the perturbation quantities. 
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IV. PERTURBATION VELOCITY 


The solution of the transformed vorticity transport equations de- 
veloped in the previous chapter can now be expressed in terms of the 
perturbation quantities P(r) and Q(r). The original functions G(r) and 
H(r) are actually components of the vector potential function W defined 


in equation (2-28), 


W=([F(r) € + G(r) €. + H(r) €,] e* (2-28) 


The axial component F({r) was previously set to zero since one function 
could be arbitrarily specified. From equation (2-27) the perturbation 


velocity vector was expressed as, 


= 


v=vuxw (2-27) 


Therefore the functions G(r) and H({r) represent the behavior of the 
perturbation velocity in the flow field. Referring back to the defini- 
tion of the curl] operation in equation (2-39), the relation for V x W 
has already been set up. The resulting expression for the perturbation 


velocity is given by, 


= = > = 
V=ue.+yve +we (4-1) 
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The components of the perturbation velocity are, 


u = dH(r) + A - 2 Gry (4-2) 
v= -ia H(r) (4-3) 
w = ia G(r) (4-4) 


Digressing for a moment, it can be readily seen that since the pertur- 
bation velocity necessarily vanishes at the pipe wall, r = 1, the 
boundary conditions at the wall in equation (3-1) can be deduced. 

Recall that the functions P(r) and Q(r) appear as the result of the 
change of variables for each specific angular wave number contained in 
subset (a) of equations (3-17) through (3-23). Analysis of the axial 
component of the perturbation velocity u becomes useful in interpreting 
the behavior of the functions P(r) and Q(r) along the radius of the 
pipe. Equations for the axial perturbation velocity in terms of P(r) 
and Q(r) are derived from equations (4-2) through (4-4) with the change 
of variable applied for each angular wave number being investigated and 


appear as follows, 


n= oO u=2Q+r 0Q (4-5) 
nea w= 3r 0 +-r20Q= 7 r P (4-6) 
Ne 60 7 eh Seon? O73 N0 =. 7 6 re -P (4-7) 


The axial perturbation velocity is computed in part VII of the main 
investigative program for each angular wave number. To prepare the 


perturbation velocity u for plotting versus pipe radius, it is first 
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normalized. Since this velocity is in general a complex quantity, the 


perturbation velocity with the largest magnitude is determined by, 


1 
ju = (u,? + ,7)° (4-8) 


cl aax 


The normalizing velocity constant will then become Ue and the maximum 


normalized velocity will become, 


= 1+ i(0) (4-9) 


a 


The remaining normalized velocities are found by dividing through by Ue. 
The normalized perturbation velocity is then plotted versus the normal- 
ized pipe radius for various axial wave numbers and Reynolds numbers at 
the three angular wave numbers investigated. The plots and results are 


discussed in detail in the chapter on Results. 
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V. NUMERICAL METHODS 


A. GENERAL METHODS USED 


In general, the vorticity transport equations are a pair of coupled 
homogeneous differential equations which are second order in the pertur- 
bation quantity G(r) and fourth order in the perturbation quantity H(r). 
A change of variables is introduced so that the number of boundary 
conditions at the axis is reduced to three for each angular wave number 
investigated, except for n = 1. Recall for this case that there are 
four boundary conditions imposed at the axis because of the introduction 
of the additional parameter H(0). By applying the change of variables 
to the original pair of equations, the vorticity transport equations are 
transformed into an equivalent pair of coupled homogeneous differential 
equations that are second order in P(r) and fourth order in Q(r). The 
solution of the vorticity transport equations will now be in terms of 
the perturbation quantities P(r) and Q(r) which are also referred to as 
eigenfunctions, and y which is an eigenvalue. The axial perturbation 
velocity is expressed in terms of P(r) and Q(r) as derived previously 
and is examined to determine the behavior of the perturbation quantities. 
The resulting eigenvalue y is of primary interest since the magnitude 
and sign of the real part of y determines the relative stability or 
instability of the pipe flow problem. 

A numerical soiution of the transformed vorticity transport equa- 
tions is now sought. The convenient matrix format of the equations is 


given below. 
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D¢p D3P D2p oP p 
[MZ] + [Mg] + [Mg] +m) 1+ Mel) § « 
n4q 029 D2q 09 | Q 


D2P oP P 
y [No] + [Nz] + [No] (5-1) 
D2Q DQ Q 


The details of the transformed coefficient matrices for n = 0, 1, and 6 
are given in Appendix B. It is convenient to approximate the deriva- 
tives of P(r) and Q(r) by a finite number of discrete unknowns along a 
radial mesh consisting of N interior points. It is not possible to 
compute the continuous values of the functions and their derivatives 
along the pipe radius using a digital computer. Therefore, discrete 
finite difference techniques are used. For the problem at hand, the 
normalized nondimensional radius of the pipe, which becomes the in- 
dependent variable, is divided into a computational mesh. The standard 
method of representing this mesh is with uniformly spaced stations over 


the interval 0 < r < 1, Figure 5-1. 


Pipe Pipe 
Axis Wal] 
r=0 =i 


ee te Ail he i 


Figure 5-1. Radial Mesh, Standard Method 
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The vorticity transport equation coefficients are then evaluated at each 
station of the radial mesh resulting in N uniformily spaced discrete 


values with N + 1 equal intervais and N + 2 total points, including 


r=0 and r=1. The spacing between the stations is denoted by 


feaa/ (N+ 1). 

A refinement of the method for generating a computational mesh is 
referred te as the half station method and is used in this numerical 
analysis. All the radial stations of Figure 5-1 are moved to mid- 


interval or one half station (h/2) toward the axis with an additional 


station appearing at one half station from the pipe wall, as shown in 


Figure 5-2. 
Pipe 
axis [V2 es 
r=0 
1 2 3 4 5 N=-4. N-3 .N-2. Nel, N 


Figure 5-2. Radial Mesh, Half Station Method 


Therefore for N interior stations, the method results in N intervals and 
N + 2 total points including the boundaries. The first station, r,, is 


a distance h/2 from the axis and the Nth station, is a distance h/2 


Py 
from the wall. The general expression for the ith station, correspond- 


ing to the ith increment of the radius, is given by, 


rae netie= 3) 304m, goals. i (5-2) 


where h = 1/N. 
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An advantage of the half station method over the standard method is 
that there is better resolution at the stations near the boundaries 
where instabilities are thought to originate and where boundary layer 
effects are a factor. Therefore, the behavior of the perturbation 
quantities can be examined more closely. Additionally, better approxi- 
mations for the derivatives of the functions P(r) and Q(r) can be 
realized with a smaller incremental distance to the points immediately 
adjacent to the boundaries. Even better resolution at the boundaries 
can be achieved by using a nonuniform computational mesh. Arnold 
{Ref. 11] introduced a change of independent variable by means of hyper- 
bolic functions which included a mesh offset parameter so that he could 
contro] the concentration of mesh points at either the axis or wall. 
The method would certainly be useful in reducing computational time by 
reducing the total number of mesh points required for a satisfactory 
solution of the governing equations. However, the full implications of 
the nonuniform mesh with respect to a and y dependence on the mesh 
offset parameter have not been fully explored. While this method is not 
used here, it may weil prove to be useful in follow-on numerical in- 
vestigations of the pipe flow problem. 

The solution for the vorticity transport equations in terms of the 
perturbation quantities P(r) and Q(r) are now represented by a finite 
number of discrete unknowns. The governing differential equations are 
rewritten in finite difference form for N stations along the mesh re- 
sulting in a finite number of unknowns namely, P,, Po, P3....P,, and Q,, 


N 


Qo, Qs. ---Qy. In general, this produces 2N simultaneous equations in 2N 
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unknowns. The details of the derivation of the finite difference equa- 
tions are explained shortly. Substituting the finite difference equa- 
tions into the vorticity transport equation matrix format of equation 
(5-1) will result in N pairs of coupled simultaneous equations. These 
equations are regrouped and subsequently soived in the following eigen- 


value problem format, 


EAT ARs = y (8) (23 (5-3) 


The column vector {X} represents the rearranged unknowns as 


= Py (5-4) 


The symbol y denotes the characteristic value or eigenvalue of the 
solution. Details of the composition of the [A] and [B] matrices are 


discussed for each specific angular wave number investigated. 
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B. FINITE DIFFERENCE EQUATIONS 


The use of finite differencing techniques is the basis for solving 
the governing vorticity transport equations in the eigensystem equation 
(5-3). It is required that the governing equations must be satisfied at 
each of the interior stations i = 1, 2, 3....N, therefore producing 2N 
equations in 2N unknowns. 

Except for stations adjacent to the boundaries, the derivatives of 
the perturbation quantities in the governing equations at the interior 
stations are usually approximated by standard central finite difference 
equations. These equations were taken from Table 7.27 of Ketter and 


Prawel [Ref. 14] and are summarized below, 


ek | 5 ie = 4ns os 
OF; = aan Ping ~ Fina * isa ~ Fis) * 35 MOF; (5-5) 
ee eee ee eee 55 N4D8F, (5-6) 
me iene } 1-2 i-1 i417 Pisa | 
er, + 1a. = ar... + BF. 5 - Ft oy nt0F, 
aece }, i=3 i Zaey thet i+1 i+2 ~ [i434 7 
(5-7) 
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where F is an arbitrary function. 
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As can be seen from the equations above, the derivative of a func- 
tion at a particular station can be approximated by a linear combination 
of the discrete unknown functions about a central point. The term 
outside the brackets in each equation represents the truncation error of 
the corresponding finite difference approximation. Notice that the 
truncation error of equations (5-5) through (5-8) are all of order h*# 
(Ch*). It was hypothesized and shown by Gawain and Ball [Ref. 15] that 
finite difference approximations of consistent order truncation error 
reduced the overall error of the numerical soiution in a typical problem 
by a factor of five for truncation error Oh?. For the present numerical 
analysis, a computational mesh of N = 50 was used for all cases. Now 
with the interval h = 1/N = 0.02, examination of the consistent fourth 
order truncation error term will yield an error of the order of 5x10” 
for all four derivatives approximated, since the magnitude of the 
normalized functions and derivatives are of order one. For an interval 
size of h = 0.02, a consistent truncation error Oh* should significantly 
reduce the overall error of the solution. 

Special problems arise in the derivation of the finite difference 
equations at stations near the boundaries. Several approaches can be 
taken to approximate the derivatives at these stations, which include 
backward, forward and central finite differencing techniques. Unfor- 
tunately, the central method previously discussed will involve imaginary 
stations outside the boundaries for stations close to the axis and wall, 
7.e., stations i = -1, -2, -3, and N+1, N+2, N*3 are involved in the 
fourth derivative approximation at station i = 1 and i = N, respectively. 


Unlike the boundary conditions for many problems well suited for this 
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type of numerical analysis, the vorticity transport equation boundary 
conditions for the pipe flow problem are very complicated by comparison. 
The boundary conditions for each angular wave number are different; some 
involve additional parameters while others contain the eigenvalue y. In 
other words, the conditions at the imaginary stations for conventional 
problems can normally be expressed in terms of conditions at the real 
stations in the computational mesh. It is therefore impractical to use 
the central method for the pipe flow problem because of the complicated 
nature of the boundary conditions. 

The method of forward differencing is used near the pipe axis and 
backward differencing is used near the wall. These two differencing 
methods will be referred to, collectively, as noncentral finite differ- 
encing techniques. It will be shown later that the two methods are 
essentially equivalent. A rigorous derivation of the noncentral finite 
difference equations is based on the Taylor series (power series) expan- 
sion of the function at a point. Since the governing equations contain 
the fourth order derivative of Q(r), the derivation for this case is 
examined in detail here. 

The Taylor series approximation of Q.(r) at station i = 1 appears as 


equation (5-9). 


Q, = Q(0) + 0Q(0) 3 + D2Q(0) S— + 03Q(0) S— + 0#Q¢0) 4 — 
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The approximation for Q: at stations i = 1, 2, 3, 4, 5, 6 are expressed 


in convenient matrix format below, 


ram 
Mle 


NI~N 


Q(0) 
hDQ(0) 
h202Q(0) 
h?D°Q(0) 
h*D4Q(0) 
h°D°Q(0) 
h®D°Q(0) 
h7D7Q(0) 


+ {E,} h®D°Q(0) 


6x8 


(5-10) 


where the dots represent matrix elements that can be easily deduced. 
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Since the fourth order derivative of Q(r) is present in the vorti- 
city transport equations, there will be two boundary conditions speci- 
fied at the axis and two at the wall. It can be shown that two columns 
of the matrix in equation (5-10) can be eliminated since two of the 
functions in the coiumn vector are either zero or can be expressed in 
terms of additional parameters. The matrix, which will be denoted as 
[AA], now becomes a square 6 X 6 matrix. Notice that the column vector 


{E,} represents the last nonvanishing term in the Taylor series. The 


approximation for the first derivative of Q;(r) at station i = 1 is 
given by, 
2 3 4 
: @) 3) @ 
DQ, = 0Q(0) + 07Q(0) 5 + D°Q(0) —— + 04Q(0) -G— + 05Q(0) 
5 6 
@) @) 
+ 0°Q(0) =e D’Q(0) BT a fees as ($=11) 


The Taylor series expansion for the second, third, and fourth deriv- 
ative can be found by successive deferentiation of equation (5-11). The 
first four derivatives of Q; at station i = 1 can be expressed in matrix 


format as, 


61 


; oo 
“tpvew eta ci gnsesag 21.00 4e ~~ 


tiny 
~iseq2 ene) tiano> yrebnved ow! o€ 


oy = 
" | f 
y 603 pucde «6 nea JL... lige altd : ed 
ey ia | | 
; : os : irs 
be. mtle at @é@> COPS) me 
7 iws epnrer te i 


- 5 


or! sara ‘© 7"o32sv veioe 
rg 5 } 


: ) >) 9398p TRO 


, & @® d @d & 
Se ae ss ee eT 
1,2 es 1.2 cS 
es) &) G) 
nO, ti. Om ee 
h2D2Q, 
( an ae 
n404Q, Dy See et Rr Pa: ee 
2 3 
oye. 
Ce © a 0 a (5) oe ap 
4x8 
Q(0) 
hDQ(0) 
h202Q(0) 
h302Q(0) 
h#D4Q(0) + {Eo} h8D#Q(0) (3-12) 
h505Q(0) 
h8D&Q(0) 
h7D7Q(0) 


It can be shown here as well that two columns of the matrix in 
equation (5-12) can be eliminated so that the matrix denoted [CC] is now 
a 4 X 6 matrix. The column vector {Ej} also represents the last non- 
vanishing term in the Taylor series. It can be seen that the approxi- 
mations for the derivatives of Q; at stations i = 2 and i = 3 can be 
determined by replacing (1/2) in the [CC] matrix with (3/2) and (5/2), 


respectively. 
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To show how the noncentral finite difference approximations near the 
axis are formulated, the case for angular wave number n 2 6 is examined. 


The boundary conditions at the axis for the function Q are, 


Q(0) = 0 and DQ(0) = 0 ($=13) 


Applying these conditions to equations (5-10) and (5-12) will eliminate 
the first two columns of the [AA] and [CC] matrices since the first two 
terms of the column vectors vanish. The truncation error terms {E,} and 
{E2} are very small quantities and are neglected at this point. In 


shorthand notation, equation (5-10) now becomes, 


Qi h2D2Q(0) 
Qo h?D3Q(0) 
Qs h*D4#Q(0) 
= [AA] (5-14) 
Q4 6x6 h°D5Q(0) 
Qs h®DQ(0) 
Qe h7D7Q(0) 


and equation (5-12) becomes, 


h202Q(0) 

hOQ,; h3D3Q(0) 

h202Q, h4D4Q(0) 
= NehCC] (5-15) 

h303Q, 4x6 h5D5Q(0) 

h404Q, h®D®Q(0) 

h7D7Q(0) 
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Equation (5-14) is multiplied through CAAJ.- and becomes, 


h2D2Q(0) G% 
h303Q(0) Qe 
h#04Q(0) a) Qs 
= [AA] (5-16) 
h5D5Q(0) 6x6 Q4 
h°D&Q(0) Qs 
h7D7Q(0) Qe 


and is then substituted into equation (5-15) giving, 


Qi 

hDQ, Qo 

h2D2Q, -1 Qs 
= GG [AA ] (5-17) 

h3030, 4x6 «G6 Qe 

h*04Q, Qs 

Qe 


Now, the fourth derivative of Q. at station i = 1 is approximated in 
terms of the discrete functions of Q at the interior mesh points with 
truncation error Oh*. The values on the bottom row of the resulting 
4X6 matrix after multiplication in equation (5-17) are the coeffi- 
cients of the column vector Q, through Q, for D4Q,. In order to have 
consistent truncation error Oh*, the last row and last column of the 
[AA] and [CC] matrices are deleted for the third derivative of Q. 
Matrix [AA] becomes a 5 X 5 and is now inverted and matrix [CC] becomes 


a3X 5. The values on the bottom row of the resulting 3 X 5 matrix are 
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the coefficients of Q, through Qs, which approximates 0°Q,. The rows and 
columns of matrices [AA] and [CC] are successively reduced to determine 
the coefficients for D?Q, in terms of Q, through Q4, and DQ, in terms of 
Q, through Q3 so that consistent truncation error Oh* is preserved 
throughout. The noncentral finite difference equation for D*Q, will now 


appear as, 
D4Q; = pa (AiQi + AaQe + AgQs + AqQy + AsQs + AcQe) + Oh# (5-18) 


where A, through Ag are the coefficients taken from the bottom row of 
the resulting matrix in equation (5-17) 

In a similar manner, by imposing the required boundary conditions, 
this process can be carried out to determine the noncentral finite 
difference equations for the first two derivatives of P(r) at several 
points near the boundaries as required. 

It is important to note here that the noncentral finite difference 
approximations near the wall can be formulated by using the same tech- 
niques as those used near the axis. Once again taking the case for 
angular wave number n 2 6, the boundary conditions at the pipe wall for 


the function Q are, 
Q(1) = 0 and DQ(1) = 0 (5-19) 


Rather than set up complete new [AA] and [CC] matrices at different 
values of r in the Taylor series expansion near the wall, assume that 


the boundary conditions in equation (5-19) are imposed at the axis, 
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r=0. Notice that the boundary conditions in equation (5-19) are now 
identical to the boundary conditions in equation (5-13). The method for 
computing the coefficients of the nencentral finite difference approxi- 
mations at the wall can be carried out as discussed previously and, in 
this case, have already been computed. t can be shown that the A. 
coefficients of equation (5-18) are applicable to 04Q,, at the wall but 
will appear in reverse order. The noncentral finite difference equation 


for 04Q,, is given by, 


1 
4 = 4 
DmQn = AF CAgQn-5 * AgQn-g * AgQy-3 + AgQy-2 + ApQy-y * AqQq2 * OF 


(5-20) 


The methed for determining the coefficients for the remaining deriv- 
atives of Q and P near the wall proceed as though the boundary condi- 
tions were imposed at the axis; the coefficients are computed, their 
order reversed and coefficients of the odd order derivatives undergo a 
sign change. Order reversal and sign changes are accomplished in part 
II of the main investigative program when the precomputed noncentral 
finite difference coefficients are read in as data. 

As can be seen from the central finite difference approximations of 
equations (5-5) through (5-8), the third and fourth derivatives of a 
function require three points either side of the central point being 
approximated. The first and second derivatives require only two points 
either side of the central point. Therefore the central finite dif- 
ference approximations are applicable only to interior stations i = 4 


through N-3 for 04Q. and D°Q., and stations i = 3 through N-2 for 07Q., 
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DQ; , 0?P., and OP... This leaves either two or three points immediately 
adjacent to the boundaries to be approximated by the noncentral finite 
difference equations just developed. 

Since the central finite difference approximations are applicable to 
the above mentioned interior mesh points and are valid for all the de- 
rivatives of P(r) and Q(r) as well, the coefficients are computed only 
once in part I of the main investigative program in double precision 
format. The coefficients of the noncentral finite difference approxi- 
mations, which are dependent on the boundary conditions for each angular 
wave number, are computed in an independent computer program. The 
resulting noncentral coefficients for the derivatives of P(r) and Q(r) 
at mesh points near the boundaries are also computed in double precision 
format and are read into part II of the main investigative program as 


data. 
C. SPECIFIC METHODS FOR n = 0, 1, AND 6 


The vorticity transport equations have now been expanded into a set 
of discrete simultaneous equations with respect to a half station 
computational mesh of N points (stations). The derivatives of the 
governing equations were approximated by central and noncentral finite 
difference coefficients, with consistent fourth order truncation error, 
in terms of the discrete unknown functions of the system. The con- 
venient matrix format of the eigensystem in equation (5-3) can now be 
solved on a digital computer using a specialized Fortran subroutine, 
which wil? be discussed later. It is important to understand the com- 


position of the [A] and [8] matrices and the column vector of unknowns 
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{X}. The problem setup is examined in detail for each angular wave 
number investigated. In the following discussions, the vorticity trans- 
port equations are referred to in the transformed state after the change 
of variables to P(r) and Q(r) have been made. The boundary conditions 
in equations (3-17), (3-18), and (3-23), and the coefficients in 
Appendix B also apply. 

For the case n = 0, the pipe flow problem becomes greatly simplified 
in that the equations become uncoupled and only the solution for Q(r) is 
investigated for reasons explained earlier. The governing equation in 
terms of Q(r) is represented by equation (2-38). At the outset, it was 
determined that a mesh size of N = 50 was sufficient to produce a satis- 
factory numerical solution of the problem. It was necessary to deter- 
mine when the solution was relatively independent of mesh size. Prelim- 
inary trial computations were performed for various mesh sizes, N = 25 
to 100. From N = 45 to 100, the solution of the equation for n = 0 did 
not vary appreciably. Additionally, the solution was relatively in- 
dependent of Reynolds numbers below 15,000 for mesh sizes of about 50. 
In order to reduce computation time and the eigensystem matrix size, a 
mesh of N = 50 was used in this analysis. 

In general for angular wave number n = 0, the [A] and [B] matrices 
are 50 X 50 and appear as banded diagonal matrices. For a consistent 
order truncation error Oh* used here, the derivative of highest order 
requires the greatest band width. From the governing equation in Q, it 
can be seen that the highest order derivative in [A] is D4Q and in [B] 
is D?7Q. At stations near the boundaries, the noncentral finite differ- 


ence coefficients require a band width of six in matrix [A] and four in 


68 


= 


a 
duce aetuons tock oe Teagan nee re 
7 - 
eead waTatrvov O42, 2nokezUaere p rRe 


4 
gis at? WeITt Bthse Soro; men? alge O: 

; Nd : 
vb T Ad al .ebam sed evert ‘70's : 


ves fe “(8reE) ‘ 


oofq ot? (Oar 
Woon! Sampo” ange 
"Tei oxS oRhvese 
sees Zi 

ict & Ja 

iso 


; ale iwi ce ont 


s'7ebAued ote 


iifsm nt <'2 90 otw td 2 er 


i 


7 


matrix [B]. At interior stations, the central finite difference co- 
efficients require a band width of seven for [A] and five for [B]. 
Therefore the [A] and [B] matrices of equation (5-3) contain nonzero 
elements as shown in Figures 5-3 and 5-4. The X elements contain non- 
central finite difference coefficients and the 0 elements contain 
central difference coefficients. The X and 8 symbols denote the station 
at which the approximation is made. The composition of the column 


vector {X} can also be seen in Figures 5-3 and 5-4. 


XXXXXX 0, 
XXXXXX 
XXXXXX Qo 
0008000 Qs 
0008000 . 
0008000 Q, 
0008000 : 
0008000 G 
XXXXXX N-2 
XXXXXX Qne1 
XXXXXX 
X Qy 


Figure 5-3. [A] Matrix for n= 0 
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XXXX Qo 
00900 Q 
00800 5 
00800 
00800 Q; 
00800 Q 
00800 ae 
00900 7 
XXXX 2 : 
KXXX N 


Figure 5-4. [B] Matrix for n =0 


The actual values of the elements of the [A] and [B] matrices are 
computed in part IV of the main investigative program. The calculations 
are rather lengthy and tedious and are accomplished in iterative loops 
in part IV. To save computation time, calculations are made only for 
the nonzero elements of the matrices. Recall] that the governing equa- 
tion in Q(r) for n = 0 was given in equation (2-43). In discrete form 


for the ith station, equation (2-43) appears as, 


Fi 4 BS 3 4 2 is 4 = 
Ma; D Q; i M3. D Q. si Mo; D Q; i Mi; DQ. i Mo; Q. 


, 2 - , 
yON2 D Q; ts Ny. DQ. + No; Q:) (5 214 


where the primed quantities are element (2,2) of the corresponding 
transformed vorticity transport equation matrices given in Appendix B, 
Part A. At station i = 1, the top row of the [A] matrix will have six 


elements and the top row of [B] will have four elements. 
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For example, the matrix elements for the first row of [A] can be 
found as follows, where the a. coefficients are analogous to the co- 
4 
efficients A. of equation (5-18) and are equal to a, = A./h . Similarly 
= 3 = 2 = 
b; B./h » C C./h , and d. 0./h, where b;, c., and d. are the 


noncentral finite difference coefficients of the remaining three 


derivatives. 
Aa 1 3} = (Ma a, + M3 Bae M2 oS ag Mi dy>* Mo) {Qi} (3°22) 
Ay 2 105} = (Ma ag + Ma. bo + Me Go t Mi do) {Qo} (5-23) 
Ay 3 {Q3} = (Ma ag + Ms bs + Mo Cs + Mi d3) {Q3} (5-24) 
Ar 4 {Q4} = (Ma aq t+ Ma D4 + Mo C4) {Q4} (5-25) 
Ars {Qs} = (Ma a5.* Mg bs) {Qs} (5-26) 
A 6 {Qe} = (Ma ag) {Q6} (5-27) 


where Ay 1 is the matrix element in the first row and the first column 
of [A], and so forth. This procedure is carried out for N equations 
along the computational mesh until all the nonzero matrix elements of 
[A] and [B] are computed. 

Derivation of the noncentral finite difference coefficients for 
other angular wave numbers is fairly straightforward when accomplished 
by the methods outlined in the previous section. The boundary condi- 


tions for Q at the axis for n = 0 are, 
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0Q(¢0) = 0 and D7Q(0) = 0 (5-28) 


After imposing these conditions, the second and fourth columns of the 
matrices in equations (5-10) and (5-12) can be eliminated since the 
corresponding terms in the column vectors vanish. The coefficients are 
then computed as described before in an independent computer program and 
then read into the main investigative program as data. Notice that the 
bondary conditions for the function Q at the wall are the same for n = 0 
and n 2 6. Therefore, the noncentral finite difference coefficients 
need to be computed only once. 

The case for angular wave number n 2 6 will be examined next since 
it more nearly represents the general case. Since the vorticity trans- 
port equations remain coupled for n 2 6 and no additional parameters are 
introduced, a system of 2N equations in 2N unknowns will result. Fora 
mesh size of 50, the [A] and [B] matrices become 100 X 100 and are 
solved in the form of equation (5-3). The first N equations of the 
system are represented by the discrete form of equation (2-37), which 
after transformation by the appropriate change of variables corresponds 
to the top equation of (5-1). The second N set of equations are repre- 
sented by the discrete form of equation (2-38), which after transfor- 
mation by the appropriate change of variables corresponds to the bottom 
equation of (5-1). As a matter of convenience for programming, the [A] 
and [3] matrices are partitioned into four 50 X 50 quadrants, which are 
individually banded diagonal matrices. The properties of the [A] and 


[B] matrices for n = 6 are summarized in Table 5-1. 
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Table 5-1. Properties of the [A] and [B] Matrices for n 2 6 


HIGHEST ORDER NONCENTRAL CENTRAL 
DERIVATIVE COEFFICIENTS COEFFICIENTS 
MATRIX | QUADRANT PRESENT REQUIRED REQUIRED 


* Number of stations near the boundaries where the noncentral finite 
difference approximations apply. 


It can be deduced from the data given in Table 5-1 that the [A] and 
[B] matrices will contain nonzero elements as shown in Figures 5-5 and 
5-6. The column vector {X} of unknowns is also given. 

The method for determining the noncentral finite difference co- 
efficients for n 2 6 was discussed in a previous example. 

The case for n = 1 is unique in that an additional parameter H(0) is 
introduced by the change of variables and a pair of coupled boundary 
conditions at the axis in equation (3-18b) contain the eigenvalue y. 
These peculiarities are significant factors in correctly setting up the 
probiem. Recall that in the coupled governing equations (2-47) the 


additional parameter H(0) appears explicitly and two additional column 
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Figure 5-6. 
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onde 


matrices are required in the system of equations. This was due to the 
change of variables introduced by equations (2-44). The appearance of 
the eigenvalue y in the pair of boundary equations at the axis comes 
about from the rigorous method in which the boundary conditions were 
derived in Chapter III. The coupled axis boundary conditions in equa- 


tion (3-18b) are repeated here, 


96 D2P(0) P(0) =< 
- 3 + 2ialCg] + jalCo] H(0) 
D2Q(0) 
(5-29) 
P( . 
=y {2 [Dg] + a? [Do] H(0) 


where the coefficient matrices are contained in Appendix C. Notice that 
the parameter H(0) appears explicitly in these two equations as well. 
Additionally, the parameters D2P(0), D2Q(0), P{<0), and Q(0) are also 
introduced. Obviously, this poses some problems since there are only 
two additional equations and five unknown parameters. H(0) can be used 
since it appears in the 2N discrete equations of the system. Q(0) is 
chosen as the second variable as a matter of convenience. The remaining 
three functions can be approximated in terms of discrete functions of P. 
in the case for D?2P(0) and P(0), and Q. and Q(0) in the case of D#Q(0). 
By referring to the matrix form of the Taylor series expansion of a 
function at a point in equation (5-10), it can be seen how the approxi- 


mations of these functions are made. For the present situation where 
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DP(0) = 0, eliminate the second column of the matrix and equation (5-10) 


will becone, 


Py P(0) 
Po h2D2P(0) 
Pa = [AA] h3D3P(0) + Oh® (5-30) 
5x5 
P, h4D4P(0) 
Pe h>D5P(0) 


Recall that the matrix in equation (5-30) is denoted [AA] as before. 


After multiplying each side of the equation by Aaa, the result becomes, 


P(O0) Pi 
h202P(0) P. 
h3D3P(0)> = [AA] 2 p,s - One (5-31) 
5x5 
h*D4P(0) P4 
hSD5P(0) Pe 


Now P(0) is approximated by the coefficients in the first row of caa] 
and D2P(0) is approximated by the coefficients in the second row of 
[aay divided by h2, both expressed in terms of P; through Ps. Notice 
that the approximation for D2P(0) has truncation error Oh* and P(Q) has 
truncation error Oh®, which is acceptable in this analysis. 

The approximation for D2Q(0) is similar except that Q(0) will appear 
explicitly in the approximations with 0Q(0) = 0. Equation (5-10) now 


becomes, 


VW 


«ew old Saag th 
pobytqiait 


/ 
h?02Q(0) 


Q: ua 

Qe. : h3D3Q(0) 
Qs 1 h4D4Q(0) 

= {Q(0)} + [AA] + Oh? 
Q4 1 6x6 h>D5Q(0) 
Qs 1 h®D&Q(0) 
Qe 1 h7D7Q(0) 
(5-32) 


After rearranging, equation (5-32) is multiplied through by Caa] 2 and 


becomes, 


h*02Q(0) Qi- Q(0) 
h?0°Q(0) Q2- Q¢0) 
h*D4Q(0) a Qs- Q(0) 
= [AA] = Ohi * (5-33) 
h°D°Q(0) 6x6 Qa- Q(0) 
h°D®Q(0) Qs~ Q(0) 
h7D7Q(0) Qe~ Q(0) 


The coefficients for D2Q(0) are taken from the first row of caa}y 2 and 
divided by h?. 02Q(0) is now expressed in terms of Q(0) and Q,; through 
Qs and will have a truncation error Oh®. 

The two additional boundary condition equations involving y can now 
be solved in the eigensystem of equations with H(Q) and Q(0) appearing 
explicitly as the two required unknowns. It is convenient to express 
the [Cg] and [Dg] matrices of equation (5-29) in slightly different 


form, 
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[Co] ry H(0) as —[Co*] {H(0)} 
and (5-34) 


[D9] rif H(0) as [Dg] {H(0)} 


where [Co*] and [Dg*] become 2 X 1 column matrices, 


3 
(z - 4i + 21a?) -j 
[C,*] = and = [Dg*] = (5-35) 
-22 + iat } 
2a2 + Re 1 


Substituting equations (5-35) into equation (5-29) yields, 


D2P(0) ) (P 0 
- 8 ter] + 2iafCs]) ‘+ ialCy*] (H(0)} 
020) (aco) 


(5-36) 


ena {2 
= y ye thes + a2 [D9*] {H(0)} 
(aco 

The solution of the eigensystem of 2N + 2 equations in 2N + 2 un- 
Knowns is now accomplished in the form of equation (5-3). The appear- 
ance of the [A] and [B] matrices are identical to those for n 2 6 in 
Figures 5-5 and 5-6 except for the addition of two columns to accommo- 
date the parameters H(0) and Q(0), and two rows to accommodate the pair 
of boundary equations in (5-36). Nonzero elements are present in the 


[A] and [B] matrices for n = 1 as shown in Figures 5-7 and 5-8. 
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Derivation of the noncentral finite difference coefficients for the 
derivatives of P and Q at the wall are straightforward as shown for the 
other angular wave numbers. The columns of the [AA] and [CC] matrices 
are eliminated based on the fact that parameters P(1), Q(1), and DQ(1) 
are all expressed in terms of H(Q), as seen in equations (3-18c). By 


imposing the boundary conditions for Q at the wall, 


Q(1) = - H(9) = and DQ(1) = 2 H(0) (5<37) 


equation (5-17) becomes, 


Qi + (1-h) H(0) 
hDQ, 2h Q> + (1-3h) H(0) 
h202Q, 0 - 3 + (1-5h) H(0) 
= {H(O)} + [CC] [AA] 
h303Q, 0 4x6 6x6 Q, + (1-7h) H(O) 
h4*D4Q, 0 Qs + (1-9h) H(0) 
Qe. +. (i-2ih) HOO) 


(5-38 


The coefficients in the resultant 4 X 6 matrix of equation (5-38) 
are in fact the same coefficients computed for n 2 6 at the wall, and 
recall] that they appear in reverse order when used at the wall. The 
details of the complete derivation of the noncentral finite difference 


coefficients near the boundaries is lengthy and not included here. 
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D. COMPUTER PROGRAM USEAGE 


The solutions for the discrete transformed vorticity transport 
equations for the three anguiar wave numbers investigated are ac- 
complished in the main investigative computer programs. Three separate 
programs were written rather than developing one program for the general 
case. The rationale for this is as follows: the changes of variables 
for each angular wave number n $ 6 are different, the boundary condi- 
tions are different, the equations uncouple for n = 0, additional para- 
meters are introduced for some, and n = 1 has a coupled boundary 
equation that contains y that must be dealt with separately in the setup 
of the problem. It would be practical to write a general program for 
angular wave numbers n 2 6, however. 

For each specific angular wave number analyzed, the main computer 


program is divided into eight parts. Their functions are summarized 


below: 

Part I The central finite difference coefficient for the derivatives 
of P and Q are computed. 

Part II The noncentral finite difference coefficients for the deriva- 


tives of P and Q at the axis and wall are read in as data. 
Problem variables, Re and a, are read in as data. 


Part III Coefficients of the vorticity transport matrix equations are 
computed along the radial mesh. 


Part IV The [A] and [B] matrix elements are computed. 


Part V The eigenvalues and respective eigenvectors are computed in 
the subroutine EIGZC and the least stable eigenvalue is 
determined. 

Part VI Verification that the least stable eigenvalue and correspond- 


ing eigenvector satisfy the governing equations. 
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Part VII The normalized axial perturbation velocity is computed from 
the least stable eigenvalue and its eigenvector. 


Part VIII The normalized axial perturbation velocity is plotted versus 
normalized pipe radius utilizing the Versatec plotter. 

The numerical analysis of the pipe flow stability problem is cen- 
tered around the subroutine EIGZC once the problem is properly set up. 
The eigensystem of equation (5-3) is solved by this subroutine which 
returns a set of N, 2N + 2, or 2N eigenvalues for n = 0, 1, or 6, 
respectively, and a corresponding set of normalized eigenvectors. Each 
eigenvalue and its eigenvector represent a solution of the governing 
equations. The primary objective of the main program and of this in- 
vestigation is to determine the least stable eigenvalue which indicates 
flow characteristics. Recall that if Yp is positive, the flow is 
unstable. For each value of a and Re chosen, the maximum algebraic 
value of Yp is found, denoted as y*, to determine if instabilities exist 
in the flow field. The corresponding eigenvector of y*, denoted as 
{X*}, is used to compute the incremental axial perturbation velocity 
along the radius. Plots of the normalized perturbation velocity versus 
normalized pipe radius have no significance in determining stability or 
instability of the flow but are indicators of the general nature and 
behavior of the flow field perturbations. Another independent program 
was written to plot Reynolds number stability contours using data taken 
from program results for n = 1 and 6. The subroutines used to make 
plots on Versatec plotter are PLOTG, standard Versatec plotting func- 


tions and TITLE1 to print titles and legends on the plots. 
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The numerical accuracy of the solution was determined by substi- 
tuting y* and its eigenvector into the governing equations. For this 
purpose, the residual error vector {¢} was computed by equation (5-39) 


and the magnitude of the error was examined. 


bed. = LAD {X*} = ys EB] KS} (5-39) 


The Fortran subroutine EIGZC is contained in the International 
Mathematical and Statistical Library (IMSL) availiable in most large 
computer system libraries utilizing Fortran compilers. Specifically, 
the numerical analysis of the vorticity transport equations in the form 
of equation (5-3) was computed on the IBM 370/3033 (assigned processor) 
system with the Fortran H Extended compiler at the U.S. Naval Post- 
graduate School. All calculations were performed in full double pre- 
cision format to improve the accuracy of the solution by minimizing 
truncation and round off errors inherently introduced by digital compu- 
tations. The main investigative programs and other auxiliary programs 
were all written in Fortran programming language to be compatible with 
all versions of Fortran compilers. The computer program listings with 


applicable data sets are included after the appendices in this paper. 
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VI. RESULTS 
A. RESULTS FOR n = 0 


Results of the pipe flow stability problem are obtained from the 
solution of the linearized vorticity transport equations by the methods 
previously described. A generalized solution of the governing differen- 
tial equations for fully developed, three-dimensional pipe flcw has been 
accomplished for three angular wave numbers, n = 0, 1, and 6. The 
characteristics of the flow field perturbations for each angular wave 
number investigated are established by the parameters Re and a. The 
purpose of this numerical investigation was to determine flow stability 
at various Reynolds numbers and axial wave numbers for each angular wave 
number. Recall that the perturbation velocity vector Vv is defined in 
equation (2-27) as the curl of the velocity vector potential W. This 
function is expressed in terms of G(r) and H(r) and a complex exponen- 
tial form e* in equations (2-28) and (2-29). After appropriate changes 
of variables were introduced, solutions were determined from the eigen- 
system of equations which consists of a set of eigenvalues and their 
corresponding normalized eigenvectors for fixed values of Re and a. The 
resulting perturbation quantities now appear in terms of the eigenfunc- 
tions P(r) and Q(r) and are represented in discrete form in the eigen- 
vector {xX} as P. and Q., where qe 1 2. Shes ss N. As can be seen from 
equation (2-30), the real part of the complex eigenvalue ¥p will deter- 


mine the growth or decay rate in time of the perturbation. Positive 
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values of Yp indicate that the flow is unstable. The eigenvalue whose 
real part has the largest algebraic value, denoted as y*, is the least 
stable and is reported in these results to show stability trends. 

Gawain [Ref. 13] reported the numerical results for angular wave 
number n = O for the solution of the governing equation in his paper. 
Those numerical results which incorporated a purely imaginary axial wave 
number and the advancements in the derivation of the boundary conditions 
at the pipe axis were duplicated in this investigation. They are 
essentially the same results except for variations in the fourth decimal 
place of the value of y*. This is due to the improved accuracy of the 
solution as a result of the numerical methods used here. The primary 
emphasis of this investigation is centered on determining the onset of 
flow instabilities at a theoretical critical Reynolds number near 1150. 
However, solutions were obtained at other Reynolds numbers to show 
stability trends and are presented in tabular and plotted form. 

Stability data for angular wave number n = O for various combin- 
ations of Re and a are given in Table 6-1. The least stable eigenvalue 
solution of the governing equation is based on a computational mesh of 
N= 50. Table entries which contain dashes indicate that the solution 
of the governing equations was not sufficiently accurate to be included 
as valid data. Two facts become very clear in examining Table 6-1. 
First of all, the flow is stable at all Re anda for n=0. Secondly, 
there is an asymptotic trend towards neutral stability for increasing 
Reynolds numbers. This second result agrees with experimental observa- 
tions in that increasing Reynolds number has a destabilizing effect on 


the flow. It can also be seen that increasing the axial wave number has 
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a stabilizing effect on the flow. Normalized axial perturbation velo- 
city plots versus pipe radius were made at various Re and a to examine 
the behavior of the perturbation quantities corresponding to the least 
stable eigenvalue solution. Figures 6-1 through 6-9, on the following 
pages, are representative of the behavior of the axial perturbation 
velocity. The activity is concentrated near the pipe axis for n = 0 and 
dies out quickly as r approaches 1.0. Additionally, the perturbation 
velocity plots were made in order to determine if mesh fineness was 
sufficient to obtain an accurate solution of the equations. Trial 
calculations were made for uniform mesh sizes up to 100 where the 


numerical results did not vary appreciably from those reported here. 
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Figure 6-1 
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Figure 6-3 
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OB. RESULTS FOR n = 1 


The numerical results for angular wave number n = 1 indicate that 
flow instabilities indeed exist. There are four factors that may 
possibly account for the flow instabilities; (1) the complicated nature 
of the coupled governing equations for n = 1, (2) the rigorous way in 
which the boundary conditions at the axis were derived and enforced, 
(3) the introduction of additional parameters into the eigensystem, and 
(4) the appearance of the eigenvalue in a pair of the boundary condition 
equations. Table 6-2 summarizes the stability data for n = 1 over a 
wide range of Reynolds numbers and axial wave numbers. Once again the 
dashed table entries indicate inaccurate results. Notice that a mesh 
size of n = 48 was used for this case. This was necessary because of 
programming constraints that required the eigensystem to be less than or 
equal to 100 X 100. It can be seen that flow instabilities exist at all 
Reynolds numbers. This occurs coincidentally with small values of qa. 
More detailed stability data is presented in Table 6-3 for axial wave 
numbers 0S a<¢1. It is clear here that the flow becomes more unstable 
for very small values of a < 1. The stability data presented in Tables 
6-2 and 6-3 is more easily interpreted as Reynolds number contours in 
the plots of GAMMA* vs. alpha. Figure 6-10 shows Reynolds number con- 
tours over a large range of alpha. It can be clearly seen that the flow 
is unstable at all Reynolds numbers for small values of axial wave 
numbers. Figure 6-11 shows a more detailed plot of the Reynolds number 
contours over a small range of a. Here it can be seen that the flow 


becomes more unstable for decreasing values of Reynolds numbers! It can 
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also be seen here that the flow becomes even more unstable for smal] 

values of a only slightly larger than zero. The reason for these 
results cannot be readily explained but it is obvious that they are in 
direct conflict with observed experimental results. When a > 1, however, 
the contours show an asymptotic trend toward neutral stability for 
increasing Reynelds numbers; in other words, increasing Reynolds number 
has a destabilizing effect on the flow which agrees with experimental 

results. A slightly different method of presenting the data can be seen 

in Figure 6-12, as alpha contours in the plot of GAMMA* vs. Reynolds 

number. This plot further emphasizes the results discussed above; the 

flow is unstable at all Reynolds numbers for small alpha, increasing 
Reynolds number has a destabilizing effect on the flow and increasing 

axial wave number has a stabilizing effect. 

Normalized perturbation velocity plots were also produced for this 
case and show some interesting trends. The behavior of the axial per- 
turbation velocity for selected values of Re and a appear in Figures 
6-13 through 6-23 on the following pages. Upon close examination of the 
perturbation velocity plots, it becomes apparent that the computational 
mesh is not sufficiently fine to adequately approximate the velocity 
function as a is increased. This problem will be discussed shortly. 

It is interesting to note that the perturbation activity shifts from 
the pipe wall to the axis for increasing values of the axial wave number 
while Reynolds number remains constant. At small values of a, the 
activity originates at the wall. At ‘larger values of a, the activity 
shifts to the axis, for the most part. At even larger values of a, the 


activity moves away from the axis to the region between the axis and 
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wall. However, some degree of activity remains at the wall for the last 
two cases. The change in the perturbation activity as described above 
can be seen in Figures 6-14 through 6-16 for Re = 575 and in Figures 
6-19 through 6-21 for Re = 2300. The shift of perturbation activity to 
the axis with increasing axial wave number deserves comment here. 
Schlichting [Ref. 16] observed that the transition to turbulence is 
"characterized by the ........ appearance of self-sustaining turbulent 
flashes with emanate from fluid iayers near the wall along the tube." 
Therefore, it appears that unstable flows are associated with pertur- 
bation activity near the wall. As the activity shifts away from the 
wall with increasing axial wave number, the flow becomes stable. Addi- 
tionally, after the activity moves to the axis, some degree of activity 
remains at the wall. While the solution appears to remain stable, it is 
recommended that this peculiarity be examined further. The use of a 


finer mesh to resolve this activity at the wall would be useful here. 
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C. RESULTS FOR n= 6 


Results for angular wave number n = 6 are similar to those for n= 0 
in that the flow is stable at all Reynolds numbers. Stability data is 
presented in Tables 6-4 and 6-5 where it can be determined that this is 
the case. The data is once again presented as Reynolds number contours 
in the plots of GAMMA* vs. alpha in Figures 6-24 and 6-25 and as alpha 
contours in the plot of GAMMA* vs. Reynolds number in Figure 6-26. The 
trends are clear; the flow is stable at all Reynolds numbers, increasing 
Reynolds number has a destabilizing effect on the flow and increasing 
the axiaj wave number has a stabilizing effect. 

Figures 6-27 through 6-36 are representative perturbation velocity 
plots for n = 6 at several Reynolds numbers and axial wave numbers. It 
can be seen in Figures 6-31 and 6-32 that the perturbation activity 


shifts from the wall to the axis for increasing a. 


118 


Q=n ot seat? of variate ets FSH ; 
| 2} pte vil Poet? evade 25 feng! Tee 4 
ta) penimeeted ad nea 7!) svonw C0 Game es | 

nt : vol os oeteszerq nlege Sane =f etab | i 

| eigte av “7? ® j 


yas » “AMMA 9S BT oft 1 


ete 2F wot® Oe 3 134 


ar. 


o 
& 2ef Yoet ava 


4 
a 7+ ™ Hy amtz {s-o 


| 
| 
| 


e £ Yovcs js 2284 


~ : yt > nt ‘ 
ai 


’ {lew @ny 


9 


"szuLod ysaw QG UO paseq aur sanjeA :az0N 


A SanpeauablLy aypqeys 3sea] Jo spuauy Aedaq uo yAMouy 


= U uaquNN aAem se_nBuy voy eqyeg ALL LGe4ys “p-9 alqey 


Tis 


oa 2a 
Pee [oe [om 
ba T908.¢- | artet-— 
: | é —- - i 


H-76 


“/* - r — = : 
4 oldet2 f2en! Yo 2bnaxT yere0 10 
¢3 


“ 


“sqyuLtod ysaw gg uo paseq aur sanjeA :ay0N 


xA San_eauaHLz ayqeys ysea] jo spuauy Aeodaq uo yyMoUy 


9 = U JaqunN aAPM ueiNBuy vos eIeG AZLLLGeIS “G-9 a1 qe] 


120 


ia eldtat2? tesel Ya 2bnesT yas! 10 


0¢ =n : | o0.0-. 


SUrE 


> i 


€S £8 0- \ F088 0- | 


: HOES 


— —- 
' 


RIVE .0O- | Mit 


—<—— 
ee _ 


Sysi..4- | 


ebG 0 | 


p2-9 aunbl4 


YHd IU 
O'S! Ge’ cl 0°O! Gol 0°S Gre 0°O 
~ _ = oe 
SlS OOE OS os 


NOISSY 4 IEUILSNN 


95><N SYNDLNOI YAIGWNN SQ IONASY YHd 1b SA *BWWH9 


*UWWHD 


2d 


in ih ay ~~ 6 


==? -yworoan: sJaaTenu 


G2-9 aunbl4 


YHd IU 
me en 04 oe!) O° 


ostt 
I 
O0Ee 


o0Sh 

a ae 

OES ao tae a ee 
NOIISY J IGbISNN 


9 = "Nh SHNOLNOI YAGWNN SQ IONASY YHd 1b SA *BWWH9 


| XHWWHD 
122 


oO 


a oe 


_ 


nvaiveare BECTON 


—_ 
a 
_ - <-> 
a % 
iS 
s 


- a 
Conia! ie 


92-9 aanbLy 


YAIGWAN SQ IONA AY 
"0002 "0009 "0005 "000h “O00E "000¢ "0001 


es 4 
sae eo 


“fh 


125 eee a JISULS 


Sel 


NOTSIY J IGULSNN 


9: =: N SYHNOLNOI YHd 1b YIGWAN SQ IONASY SA *bWWHI 


*YWWHD 


23 


NORMALIZED PERTURBATION VELOCITY 


NORMALIZED PERTURBATION VELOCITY VS RADIUS 


1.0 
N = 6 
MESH = SO 
ALPHA = QO.00 
REY NA = 575.0 
GAMMA 
3.5 
-0.0 b3o-0-3-B-0-30-0462e9ars a --9-D-D-O-B-D-D-9--1--H-9-O-D--0-O-5-0-5-5-5-5-- 2 
-0.5 
OCTAGGN = U (REAL) 
DIAMOND = U(IMAG) 
-1.0 


0.0 B.2 Q.4u 0.6 0.8 iif 
PIFE RADIUS 


Figure 6-27 


124 


es a 


surone av YTI20J3V WOT TAaRUTe 


oh, ; 


C . A 

Be a 
n an.0 = Ae 
‘ aos Ry YS? 

“AM¥RG 


NORMALIZED PERTURBATION VELOCITY 


NORMALIZED. 


l. 


Q. 


5 


= 


ad. 


.G 


PERTURBATION 


N = 
MESH = 
ALPHA = 
REY NA = 
GAMMAx = 


OCTAGON = 
DIAMOND = 


Pipe 


VELOCITY 


U (REAL) 
U (IMAG) 


Figure 6-28 


VS RADIUS 


ig 


NORMALIZED PERTURBATION VELOCITY 


NORMALIZED PERTURBATION VELOCITY VS RADIUS 


ro 
0.5 
=u..0 
=0.5 
OCTAGON = U (REAL) 
DIAMOND = U (IMAG) 
#10 


Q.0 G2 0.4 0.6 0.8 L. 
Fire RADIUS 


Figure 6-29 


126 


NORMALIZED PERTURBATION VELOCITY 


NORMALIZED PERTURBATION VELOCITY VS RADIUS 


a 


N = 6 

MESH = 50 
ALPHA = O.00 
REY NR = 1150.0 
GAMMAx 


OCTAGON = U (REAL) 
DIAMOND = U (IMAG) 


~1.0 
0.0 Wisse 0.4 0.6 0.8 l. 


BIRE BAGIUS 


Figure 6-30 


127 


Soc ome ides -_ 


i ie 
su1dRa 2v YT130/30 WBITAS 


NORMALIZED PERTURBATION VELOCITY 


NORMALIZED PERTURBATION VELOCITY VS 


1.0 
N = 6 
MESH = SO 
ALPHA = 2.00 
REY NR = 1150.0 
GAMMAx = -0.4386 
0.5 
-l. 1 eeecceascsessscesscosoos9s9504 
-0.5 
GOTAGGN = U (REAL) 
DIAMGNO = VU (IMAG) 
-1.0 


0.0 6 0.4 ite. 
rire Aeerus 


Figure 6-31 


128 


GO: 


RADIUS 


avtgaAa ey YTISOIay 


NORMALIZED PERTURBATION VELOCITY 


NORMALIZED PERTURBATION VELOCITY VS RADIUS 


£.0 
Q.5 
-0.0 
-0.5 
OCTAGON = U (REAL) 
OIAMGNO = U (IMAG) 
~1.0 


0.0 G.2 0.4 0.6 0.8 Is 
RIPE RADIUS 


Figure 6-32 


129 


orl oe er ou t 


SuloRs ey TFLIBIsV" 


NORMALIZED PERTURBATION VELOCITY 


NORMALIZED PERTURBATION VELOCITY 


1. 


0. 


) 


.2 


-1. 


N = 6 

MESH = $0 
ALPHA = Q.SO0 
REY NR = 2300.0 
GAMMAx = -0.1390 


OCTAGON = U (REAL) 
DIAMOND = U (IMAG) 


0.2 0.4 0.6 
Pare SeoTuS 


Figure 6-33 


130 


VS RADIUS 


0.8 1h 


em 2 le 


sytash ev FTIsO Fav" 


NORMALIZED PERTURBATION VELOCITY 


NORMALIZED. PERTURBATION VELOCITY VS RADIUS. 


1.0 2 
0.5 
-0.0 
-0.5 
OCTAGON = U (REAL) 
DIAMOND = U (IMAG) 
-1.0 


0.0 Gi.:2 0.4 0.6 0.8 1 ..@ 
PIPE RABIUS 


Figure 6-34 


131 


aus 


mHITaA ey YT gous 


NORMALIZED PERTURBATION VELOCITY 


NORMALIZED PERTURBATION VELOCITY VS RADIUS 


ue 


0. 


) 


=i. 


REY NA = 4600.0 
GAMMAx = -0.1712 


OCTAGON = U (REAL) 
DIAMOND = U (IMAG) 


0 G2 0.4 0.6 


FIFE RADIUS 


Figure 6-36 


132 


ee ea, 


SUTRA ZV TTIIOUaY 


NORMALIZED PERTURBATION VELOCITY 


NORMALIZED PERTURBATION VELGCITY VS RADIUS 


1.0 
N = 6 
MESH = 50 
ALPHA = 16.00 
REY NR = 4600.0 
GAMMAx = -0.7969 
0.5 
-0.0 8-B43+49-39-3-+4+3-9-3-339-3390006306306+ Bt BS SS SES 
=0.5 
GCTAGON = VU (REAL) 
DIAMOND = U(IMAG) 
| 0 1 ! i 


0.0 0.2 0.4 0.6 0.8 ib 
RIPE RADIUS 


Figure 6-36 


133 


Vikas lait i 


auToar cy iTPIOs wee 


D. NUMERICAL ACCURACY 


Numerical accuracy of the solutions of the governing vorticity 
transport equations was verified by substituting the least stable eigen- 
value and its corresponding eigenvector into the eigensystem relation 
given in equation (5-39). The error vector was examined to determine 
the residual error of the solution. For the most part the error was of 


le to 10°?, a most satisfactory result. 


the order 10° 

Upon closer examination of the error vector and the stability trends 
in the plots, it became evident that the accuracy of the numerical 
solution for the governing equations was not only dependent on the 
Reynolds number as previously mentioned but on the axial wave number as 
well. Plots of the perturbation velocity vs. radius for very large 
axial wave numbers are contained in Figures 6-37 through 6-39. The 
three figures, one for each angular wave number investigated, are 
examples of inaccurate solutions. Dependence of the solution on Rey- 
nolds number and axial wave number, became more evident as Re and a were 
increased. Rapid excursions of the perturbation velocity over a smal] 
portion of the pipe radius indicate that the computational mesh was not 
sufficiently fine to approximate the perturbation velocity function in 
the region of activity. This supposition can also be supported by 
examining the error vector where the residual error was found to be of 
the order 10 + or worse, in the region of perturbation activity. Recall 
for n = 1 that some of the activity remained at the wall as a was in- 


creased. Refer to Figures 6-20 through 6-23. The residual error was 
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nominally of the order io.> in this region. Therefore, the true repre- 
sentation of the perturbation velocity is directly dependent upon an 
accurate numerical solution of the governing equations. These inac- 
curacies can be resolved by using a nonuniform mesh in the region of 
interest as suggested by Arnold. It is suspected, however, that after a 
thorough analysis of the effects of the problem parameters and mesh size 
on the solution have been investigated, only small values of the axial 
wave number below about 10.0 and Reynolds numbers in the vicinity of 
1150 will be of primary interest. It is felt by this investigator that 
a mesh size of N = 50 and the wide range of problem parameters used here 
yielded solutions that satisfied the governing vorticity transport 


equations to a high degree and were sufficient to show stability trends. 
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NORMALIZED PERTURBATION VELOCITY VS RADIUS 
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Figure 6-37 
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NORMALIZED PERTURBATION VELOCITY 


NORMALIZED PERTURBATION VELOCITY VS RADIUS 
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Figure 6-38 
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VII. RECOMMENDATIONS AND CONCLUSIONS 


The problem of predicting the transition from laminar to turbulent 
flow of incompressible fluids in circular pipes has been pursued for 
nearly 100 years and has been presented here in a straightforward manner. 
The governing equations for the dynamics of this problem consist of the 
Navier-Stokes equation and the continuity equation. After taking the 
curl of the Navier-Stokes equation and applying the continuity principle, 
the linearized vorticity transport vector equation is obtained. The 
velocity vector potential function, which contains functions of r and 4 
complex exponential factor, was then introduced. The resulting two 
linearily independent equations form the basis for the solution of the 
problem in the form of the vorticity transport matrix equation. 

A valid solution of the governing equations was made possible be- 
cause of two factors; (1) the advancements in the linearized theory by 
introducing a purely imaginary exponential form of the axial wave number, 
and (2) the rigorous method in which the boundary conditions at the pipe 
axis were derived and those constraints enforced. After an appropriate 
change of variables was performed, the governing equations were trans- 
formed into an equivalent pair of coupled, homogeneous differential 
equations that were solved in eigenvalue problem format. Accurate 
solutions of the governing equations were achieved by the refinement and 
combination of several standard methods of numerical analysis. These 


methods include; (1) the approximation of the governing equations along 
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a computational radial mesh developed by the half-station method, (2) 
the use of central and noncentral finite difference approximations 
having consistent fourth order truncation error, and (3) the use of 
complex, double precision number representations. 

It was this author's intention from the outset to make this work as 
complete as possible. This was accomplished by including the develop- 
ment of the theory, the derivation of the boundary conditions, and a 
detailed explaination of the numerical methods used as well as obtaining 
solutions to the governing equations. This will enable future investi- 
gators to use this thesis as a starting point and as a single source 
reference for follow-on studies of this subject. 

The conclusions of this numerical analysis follow directly from 
the results of the pipe flow stability problem which were discussed 
thoroughly in the previous chapter. 

The most significant findings are: 
(1) for n = 0, the flow is stable at all Re and a. 
(2) for n = 1, the flow appears to be unstable at all Re for small 
values of a. This result is unfortunately in direct conflict 
with experimental observations. 


(3) for n= 1, unstable eigenvalues are associated with perturbation 
activity at the wall. 


(4) for n = 6, the flow is stable at all Re and a. 
The following observations are also of interest: 


(5) for n = 0, 1, and 6, increasing Reynolds number has a destabi- 
lizing effect on the flow. 


(6) for n = 0, 1, and 6, increasing axial wave number has a stabi- 
lizing effect on the flow for fixed Re. 


(7) for n = 1 and 6, increasing axial wave number causes a shift in 
the perturbation activity from the wall to the axis for fixed Re. 
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The results for angular wave numbers n = 1 and 6 are believed to be new. 
Conclusions (1), (5), (6), and (7) are in agreement with experimental 
results or previous theoretical investigations. 

It is this author's opinion that the numerical analysis of the 
vorticity transport equations is now complete for angular wave numbers 
n=0 and 6. It is felt, however, that the numerical results for 
angular wave number n = 1 warrant further study. The unresolved paradox 
is quite puzzeling in that for n = 1 the flow appears to be unstable at 
all Reynolds numbers. Continued numerical investigations for n = 1 and 
follow-on studies of other angular wave numbers should prove fruitful in 
establishing a theory which predicts the transition from laminar to 


turbulent flow of fluids in pipes. 
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APPENDIX A 
COEFFICIENTS OF THE VORTICITY TRANSPORT EQUATIONS 


The generalized coefficient matrices which appear in the basic 


vorticity transport equations (2-40) are defined below. 


0 0 
[Ma] = 3 (A-1) 
0 ih 
in 
ws 
[Mj] = + (A-2) 
3 Re 
ae 
0 ifs 
(23 os 22) = 0) 0 
1 
[Mo] = Re + (A-3) 
= i = 2a2] 0 2ia(1-r2) 
2 2 72 2 
‘s = = 2° in x + a+ =) 0 2na(2-r) 
M1] = 3 : oe 
4 in 2n2 + 3 2a2 0 pe i 
rs r r 2ia(r) 
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APPENDIX B 


COEFFICIENTS OF THE TRANSFORMED VORTICITY TRANSPORT 
EQUATIONS FOR n = 0, 1, and 6 


A. COEFFICIENTS FOR n = 0 


Since the vorticity transport equations uncouple for n = 0, only the 
solution for the eigenfunction Q(r) is sought. Therefore the details of 
the complete coefficient matrices are not required. The coefficients 
which appear here represent the matrix element (2,2) and correspond to 


the primed quantities of equation (2-43). 


1 lal rs (B-1) 
MS = - 3 (B-2) 
Ms = ( 3 + 2a?) + 2ia (r - 3) (B-3) 
My = = (3, + sa?) bie. (1 - r2) (B-4) 
4 
Mo = - a + 2ia3 (r? - r) (B-5) 
Nj =-9r (B-6) 
Ni =- 3 (B=7) 
No = a@r (B-8) 
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These coefficients also appear and are computed in part III of the main 


investigative computer program for n = 0. 


B. COEFFICIENTS for n= 1 


Since the vorticity transport equations do not uncouple for n = 1, 
the eigenfunctions P(r) and Q(r) must be solved for simultaneously. 
With the introduction of the additional parameter as a result of the 
change of variables, H(Q) appears explicitly as an unknown in the system 
of equations. This requires the addition of two coefficient matrices to 
equation (2-40). The 2 X 1 column matrix [Ms] appears on the left side 
of equation (2-47) and one 2 X 1 column matrix [N3] appears on the right 
side of equation (2-47). The complete coefficient matrices for n=1 


are defined below. 


0 0 
[M4] = r2 (B-9) 
po te 
0 ane 
[Ms] = Re (B-10) 
ae 
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1» 2 2 87 
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[Ms] = (B-11) 
- 2 ag (2022 - 19) + Zia (r2 - rf) 
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[No] = (B-15) 
0 ant 
0 ir 

[Ni] = (B-16) 
0 =r 
d+ aa? 34 

[No] = (B-17) 
rae =2 + g*r@ 
5 ia? 

[N3] = (B-18) 
a2 


These coefficients also appear and are computed in Part III of the main 


investigative computer program for n = l. 


©. COEFFICIENTS for n = 6 


The vorticity transport equations for n = 6 remain coupled. There 
are no additional parameters introduced with the change of variables to 
P(r) and Q(r). The system of equations take the form of equation (2-40) 
once again. The complete coefficient matrices for n = 6 are defined 


below. 
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These coefficients also appear and are computed in part III of the main 


investigative computer program for n = 6. 


150 


(@t-8) 


—s 
ceqge cele emnerory 


our veteqno> evil 


( 


APPENDIX C 
COEFFICIENTS OF THE BOUNDARY EQUATIONS AT THE AXIS 


The matrices which appear in the boundary equations at the axis, 


equations (3-4) through (3-8) are defined below. 


- n2 - jin 

([C,] = {C=1) 
+ 4 in - (n* + 3) 
- n2 = 2 An 

[Co] = (C22) 
+4 in = (n@ + 4) 
< n? S32 in 

2 Z 

[C3] = (C=3) 

+2in = 5 (n? + 3) 


: ) 2n? =~ (2n2 +1) 12 = dn (2 . = 
(C,] = (C-4) 


+ 4 jin (1 = = - 2 (n2 + 1) (2 = #) 


ikey§ 


2£X%a_3HT TA 2MOTT 


ait t6 enotisiee vwistedod sag ri’ "S889Q8 ~ rot 


fad benties ene (8°) noweste Gas 


+ n2 + in 
[0,4] = (C=5) 
= 2 in + (ne + 2) 


= 2 2. gh 
6 3 
[C5] = (C-6) 
21 n 
a 6 
=f ea 
[Ce] = (C=7) 
4 24 =) 
Ae = oin 
24 4 
[C7] = (C=8) 
in . e = 5) 
"64 24 
oe een 2 ayia ue _ Sia 
i“ CG a | s (3 a4 
[Cg] = (C-99 
to? an (1 = ig) = (n2 - 3) (1 - ia) 
n2 3 in 
oe Ls 
[Dg] = €C=10) 
=n + WE) 


152 


+3 
jae + 2n2 - 20°) = 2 ih 


ite1 = (C-11) 
=a an (i ooo oe > 202) 
1 0 
[Dg] = (C-12) 
0 2 
Column matrices used in equation (5-36) for n= 1. 
(a =eaa Zia?) 
(c,*] = {H(0)} (C-13) 
ee) 
(- 202 + Re) 
=A 
[D9*] = {H(0)} (C-14) 
il 


noliseuos @! bee 2957 


hen 


APPENDIX D 
SPECIAL CONDITIONS AT THE AXIS 


The determinates are formed from the corresponding matrices as 
defined in Appendix C. For angular wave numbers n > 6, all of the 


determinates become non-singular. 
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PART Vi 
VERIFY THAT THE LEAST STABLE EIGENVALUE AND CORRESPONDING 


EIGENVECTOR SATISFY THE GIVERNING PARTIAL DIFFERENTIAL 


RESIDUAL ERROR IS OUTPUT AS EPSILON. 
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PART VIEL 
PLOT THE NORMALIZED PERTURBATION VELOCITY VS PIPE RADIUS. 
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